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Abstract 

An  analysis  and  assessment  of  two  mechanisms  in  plasma  shock  interactions 
was  conducted  under  conditions  typically  encountered  in  a  weakly  ionized  glow  dis¬ 
charge.  The  mechanisms  of  a  spatially-dependent  electron  temperature  and  ad¬ 
ditional  electron  impact  ionization  at  the  shock  front  were  examined  for  effects  on 
shock  structure  and  propagation.  These  mechanisms  were  incorporated  into  an  ex¬ 
isting  one-dimensional,  time-dependent,  fluid  dynamics  code  that  uses  the  Riemann 
problem  as  a  basis  and  numerically  solves  the  Euler  equations  for  two  fluids:  the 
neutral  gas  and  the  charged  component.  The  spatial  variation  in  electron  temper¬ 
ature  was  modeled  as  a  shock-centered  rise  in  temperature.  Additional  ionization 
was  modeled  by  incorporating  a  variable  electron  temperature  and  a  quasi-kinetic 
collision  function,  for  both  unrestricted  ionization  and  ionization  mitigated  by  ion- 
electron  recombination.  Introduction  of  a  spatial  variation  in  electron  temperature 
resulted  in  a  broadening  and  strengthening  of  the  electric  field  associated  with  the 
electronic  double  layer  (EDL)  at  the  shock  front.  Results  of  unrestricted  ionization 
were  a  broadening  and  strengthening  of  the  electric  field  associated  with  the  EDL,  an 
acceleration  of  the  neutral  shock  front,  and  the  development  of  a  neutral  precursor 
ahead  of  the  shock.  Ion-electron  recombination  was  seen  to  reduce  these  effects. 
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Double  Layer  Effects 
on  Shock  Wave  Propagation 


I.  Introduction 

Interest  in  shock  waves  first  appeared  in  aeronautics  during  World  War  II. 
At  the  time,  supersonic  flight  was  thought  to  be  unobtainable  due  to  aerodynamic 
and  structural  limitations.  However,  when  dive-bombers  entered  steep  dives,  loss  of 
aerodynamic  control  was  attributed  to  the  build  up  of  shock  waves  on  aircraft  control 
surfaces  as  the  flow  over  the  wings  became  supersonic.  This  phenomenon  spurred 
serious  research  efforts  into  the  possibility  of  supersonic  flight  and  soon  after  the  war, 
the  Bell  X-l,  piloted  by  Chuck  Yeager,  traversed  the  so-called  sound  barrier.  This 
first  step  eventually  led  to  such  remarkable  developments  as  the  SR-71  and  today  s 
supersonic  fighters. 

A  new  dimension  in  supersonic  aerodynamics  may  be  emerging.  Beginning 
in  the  early  1950s,  Soviet  researchers  observed  several  phenomena  associated  with 
shock  propagation  in  weakly  ionized  gases.  Shocks  have  been  seen  to  increase  in 
velocity  ((6),  (28),  (27),  (12),  (13),  (15),  (16)),  yet  disperse  in  thickness  and  reduce  in 
strength  ((26),  (28),  (30),  (31)).  Other  structural  modifications  have  been  observed, 
such  as  an  increase  in  shock  stand  off  distance  (29)  and  the  appearance  of  a  shock 
precursor— a  region  of  elevated  density,  pressure,  and  temperature  upstream  from  the 
shock  front  ((7),  (11)).  Also,  the  aerodynamic  bodies  in  the  flow  experienced  reduced 
aerodynamic  drag  and  reduced  heating  of  their  surfaces  ((29), (22:8)).  Efforts  to 
understand  and  capitalize  on  these  effects  have  opened  a  new  field  of  study— plasma 
aerodynamics. 
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Obviously,  there  are  physical  differences  between  ionized  gases  and  unionized 
gases.  Ionized  gases  are  in  a  non-equilibrium  state,  with  neutral,  ion,  and  electron 
temperatures  differing  substantially.  Ionized  gases  also  interact  strongly  with  electric 
and  magnetic  fields.  At  the  shock  front,  electrons  diffuse  upstream  due  to  their  high 
mobility.  In  contrast,  the  ions,  being  less  mobile,  remain  more  strongly  coupled 
to  the  neutral  flow.  As  depicted  in  Figure  1,  this  charge  separation  results  in  an 
electronic  double-layer  (EDL)  and  an  electric  potential  consistent  with  Poisson  s 
equation, 


„  —  C  efn*  —  ne) 

_V20  =  V  •  E  =  -  =  — - - 

e0  e0 


(1) 


where  <j)  is  the  electric  potential,  E  is  the  electric  field  vector,  £  is  the  net  charge 
density,  e  is  the  elementary  charge,  n*  and  ne  are  the  ion  and  electron  number 
densities,  respectively,  and  e0  is  the  permittivity  of  free  space.  In  one  dimension, 
Poisson’s  equation  reduces  to 

d2<p  _  dE_  _  ejrii  -  ne) 
dx 2  dx  €0 


The  space  charge  electric  field  restricts  the  separation  of  ions  and  electrons.  The 
potential  drop,  4>,  is  related  to  the  electric  field  by  the  following: 


E  =  -V<f> 


(3) 


1.1  Problem  statement 

This  research  effort  used  analytical  and  numerical  techniques  to  characterize 
the  double  layer,  which  arises  when  a  shock  propagates  in  a  weakly  ionized  gas,  and 
assessed  double  layer  influence  on  shock  structure  and  propagation.  The  environ¬ 
ment  was  limited  to  nonmagnetized,  weakly  ionized  argon,  for  conditions  typically 
encountered  in  a  glow  discharge:  P  =  30  torr,  T  =  300  K,  Te  =  2  eV(=  23200  K),  and 
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Figure  1  Electronic  double  layer  at  a  stationary  shock  front  in  weakly  ionized  gas. 

a  fractional  ionization  of  a  =  10~6,  where  a  =  rii/nn.  Previous  treatments  simplify 
the  problem  by  omitting  several  physical  mechanisms.  Some  of  these  restrictions 
were  lifted  in  order  to  investigate  the  contributions  of  these  mechanisms  upon  shock 
structure  and  propagation.  Analytical  and  numerical  techniques  were  employed  in 
a  self-consistent  solution  of  the  resulting  Euler  equations.  Specific  attention  was 
given  to  the  contribution  of  the  electric  field  generated  by  charge  separation  at  the 
shock  front.  The  self-consistent  calculation  of  the  field  was  also  attempted  in  this 
research;  the  weakly  ionized  gas  was  represented  as  a  set  of  three  fluids  coupled  by 
collisions  and  the  electric  field.  This  would  have  employed  analytical  and  numerical 
techniques  in  a  self-consistent  solution  of  the  resulting  Euler  and  Poisson  equations. 
Details  of  this  development  can  be  found  in  Appendix  B. 


1.2  Importance  of  research 

Being  a  relatively  new  field  of  study,  the  benefits  of  plasma  aerodynamics  have 
yet  to  be  fully  tapped.  Hilbun  (22:17)  mentions  several  possible  applications: 

•  Drag  and  heat  transfer  reduction  of  hypersonic  vehicles:  Several  researchers 
have  reported  these  effects  ((8),  (29)).  Figure  2  depicts  an  advanced  hypersonic 
vehicle  concept  that  would  employ  the  benefits  of  plasma  aerodynamics  (21). 
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Figure  2  AJAX  hypersonic  vehicle  concept  (Ref:  (21)). 

•  Aircraft  without  control  surfaces:  Control  moments  are  generated  by  applying 
magnetic  torques  to  the  ionized  fluid. 

•  Boundary-layer  control:  Laminar  flow  could  be  maintained  for  conditions  in 
which  turbulent  flow  or  separation  is  expected  by  applying  magnetic  torques 
to  the  ionized  flow. 

•  Reduction  of  radar  cross  section:  For  average  flight  conditions,  radar  signals, 
which  are  on  the  order  of  10  GHz,  could  possibly  be  attenuated  by  maintaining 
a  weakly  ionized  environment  around  an  aircraft  with  fractional  ionizations  as 
low  as  10-8. 

•  Maintaining  RF  communication  with  spacecraft  during  reentry:  Communica¬ 
tions  blackout  is  typically  experienced  due  to  high  plasma  frequencies  around 
the  spacecraft.  The  presence  of  a  magnetic  field  yields  different  electromag¬ 
netic  wave  modes  for  a  given  frequency  (37).  Magnetic  control  of  the  reentry 
plasma  could  allow  for  non-attenuated  wave  modes. 

1.3  Scope  and  limitations 

The  shocks  investigated  are  one-dimensional,  ideal  processes.  Shock  phenom¬ 
ena  are  multi-dimensional  problems;  however,  limiting  the  problems  to  one  spatial 
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dimension  allows  normal  shock  relations  to  be  used,  which  will  be  discussed  in  Chap¬ 
ter  II.  At  the  shock  front,  numerous  viscous  mechanisms  occur.  Since  the  research 
objective  focuses  on  the  motion  of  charged  particles  in  the  general  area  of  the  shock, 
these  viscous  effects  are  ignored.  Therefore,  the  ideal  equation  of  state  can  be  used: 

P  =  p—T  (4) 

m 

where  P  is  the  pressure,  p  is  the  density,  ks  is  Boltzmann’s  constant,  m  is  the 
molecular  mass  of  the  gas,  and  T  is  the  temperature.  The  ideal  assumption  is  valid 
since  it  has  been  experimentally  determined  that  gases  at  low  pressures  (760  Torr  or 
less)  and  high  temperatures  (273  K  or  more)  behave  in  ideal  fashion  (4:15).  Some 
researchers  have  suggested  that  perhaps  weakly  ionized  gases  behave  in  a  non-ideal 
fashion  ((32)  and  (35)).  Saeks  and  Kunhardt  (35)  derive  an  equation  of  state  that 
incorporates  electrostatic  pressure  terms  due  to  the  plasma  components.  From  this 
equation  of  state,  an  equation  of  the  sound  speed  can  be  derived.  Mishin  (32) 
proposes  a  similar  equation  of  state,  yet  his  non-ideal  terms  are  based  on  empirical 
calculations  rather  than  a  derivation  employing  Maxwell’s  equations. 

Although  much  of  the  current  research  is  conducted  in  glow  discharge  tubes, 
this  research  effort  will  focus  on  shock  tube  problems,  also  known  as  Riemann  prob¬ 
lems,  as  the  instrument  of  investigation.  A  review  of  shock  tube  dynamics  will 
occur  in  the  next  chapter.  In  a  glow  discharge  tube,  an  electric  field  is  applied  to 
the  gaseous  medium,  which  causes  electronic  excitation  (hence  the  glow)  and  ion¬ 
ization.  Ion  density  is  maintained  via  diffusional  losses  to  the  tube  walls.  In  the 
shock  tube,  no  initial  electric  field  is  present;  however,  it  will  be  assumed  that  the 
fractional  ionization  is  constant  as  well. 

In  addition  to  the  presence  of  an  ambient  electric  field,  shock  structure  differs 
between  those  of  shock  tubes  and  discharge  tubes.  The  shocks  in  discharge  tubes 
are  generated  by  a  spark  detonation-a  momentary  burst  of  energy.  In  shock  tubes, 
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they  axe  generated  by  a  large  pressure  reservoir  and  are  consistent  with  aerodynamic 
shocks;  the  drive  of  the  aerodynamic  vehicle  provides  the  pressure  reservoir  similar 
to  that  of  a  shock  tube.  This  begs  the  question  as  to  whether  different  apparati  of 
investigation  yield  different  results.  For  a  wide  variety  of  gaseous  media,  Gorshkov, 
et  al.,  (19)  used  a  variety  of  plasmas  to  measure  shock  wave  propagation-the  plas¬ 
mas  of  glow  discharges,  pulsed  discharges,  and  decaying  glow  discharges-and  found 
that  “regardless  of  the  kind  of  gas  and  the  type  of  discharge,  the  behavior  of  the 
electron  component  in  the  region  of  the  shock  wave  front  had  the  same  characteristic 
features”.  In  addition,  similar  results  were  obtained  by  researchers  ((8),  (29))  while 
conducting  tests  in  a  completely  different  experimental  setup-aerodynamic  bodies 
moving  through  a  weakly  ionized  media. 

In  real  gases,  there  is  a  plethora  of  energy  sinks  to  consider  in  simulated  hydro- 
dynamic  phenomena.  Some  of  these  sinks  are  molecular,  such  as  rotation,  vibration, 
and  dissociation.  Other  sinks  are  both  molecular  and  atomic,  such  as  electronic  ex¬ 
citation  and  ionization.  In  the  case  of  the  present  research,  the  molecular  energy 
sinks  were  eliminated  from  consideration  by  studying  a  monatomic  gas:  argon.  To 
further  simplify  the  problem,  electronic  excitation  was  ignored-argon  atoms  are  ei¬ 
ther  ionized  or  not,  and  energy  gained  is  transferred  to  either  ionization  or  kinetic 
energy. 

1.4  Research  Approach 

Chapter  II  briefly  covers  the  background  material  of  plasma  aerodynamics. 
Before  the  effects  of  plasma  can  be  studied,  a  brief  review  of  the  basics  of  one¬ 
dimensional,  compressible  fluid  dynamics  in  the  absence  of  plasma  is  provided. 
Shock  tubes  are  also  afforded  a  brief  description.  Then,  the  phenomena  associated 
with  shock  wave  propagation  in  weakly  ionized  gases  are  identified  and  described. 
Finally,  possible  physical  mechanisms  behind  these  phenomena  are  identified:  vi- 
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brational  relaxation,  thermal  inhomogeneities,  ion-acoustic  wave  damping,  and  ad¬ 
ditional  ionization  at  the  shock  front. 

Some  of  the  previous  analytical  treatments  of  shock  propagation  in  weakly 
ionized  gases  are  illustrated  in  Chapter  III.  The  works  of  Avramenko,  et  al.,  (5) 
and  Hilbun  (22)  serve  as  the  main  fare.  As  with  any  analytic  treatment  of  complex 
physical  phenomena,  these  works  incorporate  a  number  of  restrictions  that  may 
preclude  an  accurate  evaluation  of  the  electric  field  and  its  affect  on  shock  structure 
and  propagation.  These  restrictions  are: 

•  Two-fluid  approximation,  The  charged  components  of  the  plasma  are  com¬ 
bined  into  one  fluid  by  assuming  that  the  electron  momentum  is  steady-state 
for  laboratory  conditions.  This  restriction  prevents  ions  and  electrons  from 
moving  independently  of  each  other,  which  is  necessary  in  order  to  fully  charac¬ 
terize  the  double  layer  influence  on  shock  structure  and  propagation  in  weakly 
ionized  gases. 

•  Approximation  of  electric  field.  This  is  a  by-product  of  the  two-fluid  approx¬ 
imation.  With  steady-state  electrons,  the  electric  field  balances  the  pressure 
gradient  force  and  momentum  coupling,  which  allows  the  electric  field  to  be  ap¬ 
proximated  by  means  of  the  electron  pressure  gradient,  as  opposed  to  Poisson’s 
equation  (Equation  2). 

•  Constant  electron  temperature.  This  is  assumed  to  be  the  case  in  glow  dis¬ 
charges  due  to  the  extremely  short  relaxation  time  of  the  electrons,  meaning 
any  energy  gained  is  quickly  lost. 

•  Static  neutral  profiles.  In  the  work  of  Avramenko,  the  physical  properties 
of  the  neutrals  near  the  shock  front  are  specified,  which  prevents  the  charged 
components  from  affecting  the  neutrals. 

•  Constant  fractional  ionization.  There  are  no  production  or  loss  terms;  the 
plasma  cannot  decay  nor  can  additional  ionization  occur  at  the  shock  front. 
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•  Ideal  behavior.  Electrostatic  pressures  are  prevented  from  affecting  the  gas 

pressure. 

Two  of  these  restrictions  are  lifted  in  Chapter  IV.  Research  indicates  that 
electrons  undergo  heating  at  the  neutral  shock  front  as  they  pass  through  the  EDL, 
therefore,  variable  electron  temperature  at  the  shock  front  is  investigated.  Electron 
energies  for  equilibrium  laboratory  conditions  are  generally  lower  than  the  ioniza¬ 
tion  potentials  of  the  gaseous  media  studied.  Research  indicates,  however,  that 
additional  ionization  may  occur  at  the  shock  front;  therefore,  additional  ionization 
is  investigated,  as  well.  For  both  cases,  the  numerical  results  from  modifications 
to  Hilbun’s  plasma  code  are  examined  for  effects  on  the  EDL,  shock  structure,  and 
shock  propagation.  This  research  effort  also  attempted  to  lift  the  two-fluid  ap¬ 
proximation;  however,  the  difference  in  computational  time  step  for  electron  motion 
compared  to  heavy  particle  motion  proved  to  be  prohibitively  small  to  be  accom¬ 
plished  within  the  time  and  computational  limitations  imposed  upon  the  present 
research.  The  development  of  this  attempt  can  be  found  in  Appendix  B. 

1.5  Product  of  research 

This  research  has  produced  a  modification  to  Hilbun’s  one-dimensional  plasma 
code  (22:27)  that  makes  allowances  for  variable  electron  temperature  and  additional 
ionization  at  the  shock  front.  The  aforementioned  restrictions-constant  electron 
temperature  and  constant  fractional  ionization-are  lifted.  Numerical  data  from  this 
code  modification  can  provide  a  non-restricted  data  set  to  compare  to  the  restricted 
analytical  and  numerical  solutions  of  Avramenko  and  Hilbun.  This  comparison  can 
illuminate  any  significant  contributions  to  shock  structure  and  propagation  from 
the  double  layer  electric  field  as  enhanced  by  variable  electron  temperature  and 
additional  ionization  over  a  range  of  discharge  parameters  common  to  plasma  aero¬ 
dynamics. 
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II.  Background 


2. 1  Compressible  flow  in  neutral  gas 

Compressible  flow  is  that  branch  of  fluid  dynamics  that  concerns  itself  with 
fluid  media  in  which  the  density  can  no  longer  be  considered  constant.  An  important 
phenomena  that  occurs  due  to  compressibility  is  the  shock  wave.  A  full  treatment 
of  the  shock  wave  dynamics  is  beyond  the  scope  of  the  present  research.  However, 
Hilbun  stated  it  very  simply:  “Shock  formation  can  be  understood  by  considering 
the  nonlinear  terms  in  the  fluid  equations.  These  nonlinear  terms  cause  portions  of 
the  wave  with  a  larger  amplitude  to  travel  at  a  higher  velocity  than  portions  of  the 
wave  with  a  smaller  amplitude.  Thus  the  wave  front  steepens  until  it  becomes  multi¬ 
valued  and  ‘breaks’  (22:4)”.  The  width  of  a  typical  shock  front  for  aerodynamic 
conditions  is  on  the  order  of  10-4  cm,  or  several  mean  free  paths. 

2.1.1  One  dimensional  flow.  Shock  waves  in  any  environment  are  three 
dimensional  structures.  Although,  in  some  cases,  it  is  feasible  to  model  them  in  one 
dimension.  One  such  case  would  be  at  the  nose  of  a  supersonic  aircraft.  When 
the  shock  is  modeled  in  one  dimension,  the  shock  is  normal  to  the  free  stream  and 
is  called  a  normal  shock.  Motion  in  wind  tunnels,  shock  tubes,  and  gas  discharge 
tubes  generally  assume  constant  area  flow  and  the  normal  shock  relations  outlined 
below  can  be  applied  to  them. 

Most  compressible  fluid  texts  provide  a  development  of  normal  shock  relations. 
A  brief  review  is  warranted  here  with  the  help  of  Anderson  (4) .  In  order  to  investi¬ 
gate  any  fluid  flow  process,  it  is  necessary  to  begin  with  the  continuity,  momentum, 
and  energy  equations.  In  the  analysis  of  normal  shocks,  the  flow  is  assumed  to  be 
steady  and  inviscid,  meaning  that  there  are  no  losses  due  to  friction,  diffusion,  and 
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thermal  conduction.  The  continuity  equation  is  then 


d  d 

dtp+uIxp  =  0 


Pl«l 


P2U2 


(5) 


where  p  is  the  fluid  density,  u  is  the  free  stream  fluid  velocity,  and  the  subscripts  1 
and  2  represent  conditions  upstream  and  downstream  from  the  shock,  respectively. 
This  simply  states  that  the  mass  flux  on  either  side  of  the  shock  is  constant.  As  in 
the  case  with  normal  shocks,  when  u\  >  u2,  then  p2  >  P\-  The  momentum  equation 
is  given  by 

I  W  +  4h  +  Sp  =  0 

Pi  +  Piu\  =  P2  +  P2u\  (6) 


where  P  is  the  fluid  pressure.  The  quantity  pu 2  is  often  referred  to  as  the  dynamic 
pressure.  The  energy  equation  is  given  by 


+-£(*•)+£ 

?±  +  *l  +  wi  + 

2  pi 


0 


Uo 


-±  +  —  +V)2 

2  P2 


(7) 


where  Q  is  the  heat  addition  and  w  is  the  thermal  energy,  given  by 


W  =  CyT  = 


1  keT 
7  —  1  m 


(8) 


where  7  is  the  ratio  of  specific  heat  at  constant  volume  to  the  specific  heat  at  constant 
pressure,  cy/cp,  which  is  never  more  than  5/3.  Incorporating  the  definition  of  the 
thermal  energy  and  the  equation  of  state,  Equation  4,  the  energy  equation  becomes 


^+_x_m+q=!+-^ 

2  7  —  1  m  2  7  —  1  m 
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The  speed  of  sound  is  an  important  parameter  in  the  analysis  of  compressible 
flow.  Sound  waves  impart  small  deviations  in  the  fluid  parameters;  therefore,  the 
propagation  of  a  sound  wave  is  an  isentropic  process.  A  perturbation  analysis  of 
Equations  5  and  6  yields  the  following  isentropic  relationship  for  the  speed  of  sound: 


a  = 


(10) 


where  s  indicates  an  isentropic  process.  Applying  the  ideal  equation  of  state,  Equa¬ 
tion  4,  and  the  following  isentropic  relation, 

El  =  (Ely  (li) 

Pi  Pi 


Equation  10  becomes 


(12) 


Note  that  the  sound  speed,  depending  on  the  value  of  7,  is  about  two-thirds  to  three- 
quarters  of  the  average  molecular  velocity  from  kinetic  theory,  ^ ave  ^8kBT/7rm. 

Another  important  parameter  related  to  the  sound  speed  is  the  Mach  number. 


M  =  u/a  (13) 

The  Mach  number  is  also  a  measure  of  the  directed  energy  of  the  fluid  flow  compared 
to  the  random  molecular  motion,  as  shown  by  the  following  ratio  of  the  directed 
energy  to  the  thermal  energy: 

u2/ 2  u2/ 2  _  u2/ 2  _  WJ2  _  7(7-l)  ,f 2  fl4) 

oVO 7-1)  2 
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Obviously,  M  =  0  means  that  there  is  no  directed  energy,  as  M  approaches  oo, 
almost  all  the  fluid’s  energy  is  directed,  and  M  =  1  means  that  the  flow  velocity 
exceeds  the  thermal  velocity  by  a  factor  of  -/y. 

From  these  simple  relationships  come  the  normal  shock  relations,  which  define 
flow  variable  ratios  across  the  shock  front.  For  supersonic  flow,  the  Prandtl  relation 
states  that  Mi  >  1  and  M2  <  1.  From  the  Prandtl  relation  and  the  energy  equation, 
the  post-shock  Mach  number  is  given  by 


1  +  ^M* 


7  Ml  -  ^ 


(15) 


The  following  ratios  can  also  be  derived  from  the  energy  equation: 


27 

7+1 


{Ml  -  1) 


(16) 


P2  _  «i  _  (7  +  1)M? 

Pi  ~  u2  ~  2  +  (7  -  1  )Ml 


(17) 


Zk  =  P2P1  (18) 

Ti  Pi  p2 

Note  that  all  post-shock  flow  variables  can  be  determined  by  the  upstream  flow 
variables  and  Mach  number. 

2.1.2  Shock  tube  relations.  An  important  instrument  for  studying  normal 
shocks  is  the  shock  tube.  The  shock  tube  is  a  tube  of  constant  cross  sectional  area 
closed  at  both  ends,  with  a  high  pressure  reservoir  (Region  4)  separated  from  a  low 
pressure  region  (Region  1)  by  a  diaphragm,  as  illustrated  in  Figure  3.  Regions  4  and 
1  are  also  referred  to  as  the  driver  section  and  the  driven  section,  respectively.  The 
gases  in  these  regions  can  also  have  different  temperatures  and  molecular  masses. 
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When  the  diaphragm  is  broken,  a  shock  propagates  into  Region  1  with  speed  c,  and 
expansion  waves  propagate  into  Region  4,  lowering  Pa  to  P3,  as  shown  in  Figure  4. 
As  the  shock  plows  into  Region  1,  it  increases  the  pressure  from  Pi  to  P2,  and  induces 
mass  motion,  u,  in  Region  2.  The  interface  between  the  driver  and  driven  section 
is  called  the  contact  discontinuity,  which  also  propagates  with  velocity  u.  The  flow 
field  in  the  shock  tube  after  the  diaphragm  is  broken  is  completely  determined  by 
the  initial  conditions  in  Regions  1  and  4.  The  following  equation  relates  the  initial, 
quiescent  conditions  to  the  dynamic  conditions: 

Pa  _  P2  q  _  (74  ~  _ ^-274/(74-!)  (19) 

Pl  ~~  Pl  ^/27i(27i  +  (71  +  1)(§  -  ^ 

which  gives  the  pressure  ratio  across  the  shock  front,  P2/P1,  as  an  implicit  function 
of  the  initial  pressure  ratio  across  the  shock  tube  diaphragm,  Pa/ Pi  •  After  a  period 
of  time  has  passed  since  the  initial  onset  of  the  shock,  the  shock  propagation  in  the 
Riemann  problem  becomes  relatively  constant;  therefore,  the  normal  shock  relations 
can  be  used  and  the  resulting  Mach  number  is  easily  determined  from  Equation  16. 
In  spite  of  the  interesting  effects  that  occur  in  the  shock  tube,  the  present  research 
focuses  only  on  the  shock. 

2.2  Observations  in  Weakly  Ionized  Gases 

When  shocks  propagate  in  weakly  ionized  gases,  a  number  of  phenomena  have 
been  observed  that  deviate  from  what  is  observed  in  neutral  gas  flows.  Four  of  these 
phenomena  are  investigated  in  this  work: 

•  Shock  waves  in  plasma  have  an  anomalously  high  propagation  velocity.  In 
some  cases  when  the  shock  entered  the  plasma,  the  shock  velocities  nearly 
doubled  ((6),  (28),  (27),  (12),  (13),  (15),  (16)). 
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Figure  3  Initial  shock  tube  conditions  (Ref:  (4)). 
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Figure  4  Shock  tube  dynamics.  The  shock  moves  into  the  quiescent  gas  in  Region 
1  at  speed  c.  The  contact  discontinuity  moves  with  the  fluid  behind  the 
shock  at  speed  u.  The  expansion  waves  propagate  into  Region  4  (Ref: 

(4))- 
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•  Significant  broadening  and  dispersion  of  the  shock  front  occurs  in  weakly  ion¬ 
ized  gases.  The  steep  gradients  in  the  fluid  flow  variables  that  define  the  shock 
are  significantly  reduced  ((26),  (30)). 

•  A  precursor  exists  ahead  of  the  shock  wave.  Precursors  are  regions  of  elevated 
pressure  that  precede  the  shock  by  millimeters  or  centimeters  (  (7),  (19),  (11)). 

•  Shock  strength  is  reduced.  In  some  cases  the  cross  shock  pressure  ratio  was 
cut  in  half  ((28),  (31)). 

In  addition  to  these  basic  aspects  of  shock  modification  in  a  weakly  ionized  gas, 
there  are  three  phenomena  that  are  not  examined  due  to  the  nature  of  the  model 
used  in  this  research: 

•  Modification  of  aerodynamic  drag  ((8),  (18)). 

•  Reduced  heat  flux  to  aerodynamic  surfaces  (22:8). 

•  Increase  in  shock  stand-off  distance  (29). 

The  possible  physical  mechanisms  behind  these  observed  phenomena  are  de¬ 
scribed  below.  They  are  post-shock  vibrational  relaxation,  thermal  inhomogeneities 
in  the  plasma,  and  ion-acoustic  wave  damping.  As  a  shock  propagates  into  a  region 
of  the  gas  that  is  vibrationally  excited,  the  gas  gives  up  vibrational  energy  to  the 
flow,  which  accelerates  the  shock.  Post-shock  vibrational  relaxation  is  obviously  not 
a  contributor  to  the  observed  phenomena  in  monatomic  gases,  such  as  argon. 

2.2.1  Thermal  inhomogeneities.  The  phenomena  mentioned  above  are  also 
observed  in  monatomic  gases.  Researchers  have  suggested  that  thermal  inhomo¬ 
geneities  in  the  weakly  ionized  gas  to  be  the  chief  cause  of  the  effects  ((2),  (14)). 
Elevated  temperatures  can  be  the  by-product  of  plasma  generation  or  decay.  It  is 
known  that  shocks  move  faster  in  media  with  higher  temperatures,  as  illustrated  in 
Figure  5,  where  Tx  and  Vi  represent  the  initial  gas  temperature  and  shock  velocity, 
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Figure  5  Acceleration  of  a  Mach  2  shock  as  a  function  of  thermal  inhomogeneity 
(Ref:  (22)). 


respectively,  and  T2  and  V2  represent  the  gas  temperature  and  shock  velocity,  re 
spectively,  in  a  region  where  T2  >  T\*  The  increase  in  shock  velocity  is  primarily  due 
to  the  fact  that  the  sound  speed  increases  as  temperature  increases  (Equation  12). 
The  shock  accelerates  because  it  encounters  an  area  in  which  the  gas  molecules  are 
of  higher  energy;  therefore,  less  energy  from  the  shock  is  needed  to  accelerate  the 
upstream  molecules  and  the  shock  wave  accelerates.  Hilbun  devoted  a  large  portion 
of  his  research  to  thermal  inhomogeneities  and  developed  a  two-dimensional,  time- 
dependent  thermal  code  to  investigate  this  mechanism  (22:28).  The  thermal  cause 
has  its  detractors  though  ((26),  (28),  (6),  (27)),  who  have  shown  that  shock  velocities 
in  plasma  exceed  those  that  would  be  expected  from  a  solely  thermal  configuration. 


2.2.2  Ion-acoustic  wave  damping.  According  to  Jones  (24:77),  ion-acoustic 
waves  are  purely  electrostatic  waves  in  an  nonmagnetized,  ion-electron  plasma.  The 
dispersion  relation  for  such  waves  is  given  by 
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where  u  and  k  are  the  frequency  and  wave  number  of  the  disturbance,  respectively, 
ue  and  Ui  are  the  electron  and  ion  plasma  frequencies,  respectively,  and  ve  and  v, 
are  the  electron  and  ion  thermal  velocities,  respectively.  It  is  assumed  that  a 
and  that  Te  >  Th  which  leads  to  vt  <  u/k  <C  ve,  so  that  the  dispersion  relation 

reduces  to 


1  + 


U>7 


u2  -  k2v2 


=  0 


(21) 


If  it  is  also  assumed  that  u>  cjj,  and  after  some  algebraic  manipulation,  the  group 
velocity  can  be  extracted  from  the  dispersion  relation: 


(^)2  =  -uf  +  ^(— )2  (22) 

vfc'  1  V 

The  plasma  frequency  is  given  by  u2  =  nPq2/eomp,  where  np  is  the  particle  number 
density,  and  q  is  the  charge  on  the  particle.  After  applying  Equation  12  to  both  ve 
and  Vi,  and  if  it  is  also  assumed  that  the  plasma  is  quasi-neutral  and  singly  ionized, 
then  the  ion-acoustic  velocity  is  given  by 


LU 

k 


—  Via 


(23) 


where  7e  is  taken  to  be  unity.  According  to  Jones  (24:78),  the  ion-acoustic  wave  is  a 
compression  wave  in  a  plasma  analogous  to  an  ordinary  sound  wave  in  air,  where  the 
ions  provide  most  of  the  inertia  of  the  wave,  while  the  electrons  provide  the  pressure 
to  drive  the  wave.  Note,  however,  that  via  can  be  an  order  of  magnitude  larger  than 
a,  depending  upon  the  ratio  Te/Tf. 

Several  researchers  ((26),  (6),  (7),  (28))  have  suggested  that  it  is  the  damping  of 
these  waves  that  accelerate  shocks  in  plasma.  Mishin  (28)  stated  it  most  eloquently. 
“The  high  effective  sound  velocity  in  the  plasma-several  times  the  thermal  sound 
velocity-can  be  explained  by  assuming  that  an  intense  mechanism  is  operating  to 


17 


convert  some  of  the  energy  of  the  shock  wave  into  kinetic  energy  of  the  neutral 
particles  ahead  of  the  shock  front.”  In  the  thermal  case,  the  shock  accelerates 
because  it  encounters  an  area  in  which  the  molecules  are  of  a  higher  energy ,  therefore, 
less  energy  from  the  shock  is  needed  to  move  the  molecules  and  the  shock  accelerates. 
In  the  case  of  plasmas,  the  electric  field  with  ion- neutral  momentum  coupling  serves 
to  energize  the  particles  ahead  of  the  shock,  and  the  shock  does  not  expend  as 
much  energy  to  move  through  the  gas.  Basargin  and  Mishin  (6)  state  that  the 
shock  propagation  velocity  is  a  hybrid  of  the  non-plasma  shock  velocity  and  the 
ion-acoustic  wave  velocity. 

2.2.3  Additional  ionization.  Some  researchers  indicate  that  additional 
ionization  may  contribute  to  the  observed  plasma  effects  (19).  As  electrons  are 
accelerated  by  the  space  charge  field,  a  small  fraction  of  them  reach  kinetic  energies 
on  par  with  the  first  ionization  potential  of  the  gas.  The  charged  particles  produced 
at  the  shock  front  are  generally  restricted  to  that  location  due  to  the  local  fluid 
dynamics,  the  electric  field,  and  the  mutual  attraction  of  the  charged  particles.  The 
net  effect  of  this  ionization  is  an  exponential  increase  in  charged  particles  with  time, 
which  increases  the  local  fractional  ionization  at  the  shock  front,  making  ion- acoustic 
wave  damping  more  effective.  This  mechanism  is  explored  in  depth  in  Chapter  IV. 


18 


III.  Literature  Review 


The  effects  of  plasma  on  shock  wave  propagation  have  been  studied  since  the  1950s. 
In  1965,  Jaffrin  (23:616)  attempted  to  numerically  solve  the  Navier-Stokes  equations 
for  a  ternary  fluid  of  neutrals,  ions,  and  electrons.  Due  to  the  complexity  of  the 
equations  and  the  limited  computational  resources  at  that  time,  all  three  fluids  were 
combined  into  one  in  order  to  facilitate  a  solution.  His  solution  for  weakly  ionized 
gases  yielded  a  charged  precursor  upstream  from  the  shock,  but  the  neutral  flow  was 
unaffected.  The  analytical  foundation  of  shock  wave  structure  in  weakly  ionized 
gases  was  developed  in  the  early  1980s  by  Avramenko,  Rukhadze,  and  Teselkin  (5) 
and  is  therefore  reviewed  in  detail.  Hilbun  extended  their  steady-state  approach 
(22:69)  and  continued  with  a  three-pronged,  time-dependent  approach,  examining 
plasma  effects,  post-shock  vibrational  relaxation,  and  thermal  inhomogeneities.  The 
present  research  is  an  extension  of  Hilbun’s  investigation  of  plasma  effects. 

3.1  Steady  State  Treatment 

3.1.1  Analytical.  Avramenko,  et  al.,  (5)  sought  an  analytical  solution  to 
the  problem  of  shock  propagation  in  weakly  ionized  gases,  particularly  focusing  on 
the  flow  variables  of  the  charged  components  thereof.  In  order  to  accomplish  this 
feat,  they  applied  all  of  the  restrictions  described  in  Section  1.4.  These  restrictions 
were  the  two-fluid  approximation,  the  electron  pressure  gradient  approximation  of 
the  electric  field,  constant  electron  temperature,  static  neutral  profiles,  and  no  pro¬ 
duction  or  loss  terms.  The  treatment  further  imposes  the  following  restrictions: 

•  The  ion  temperature  is  constant  through  the  shock. 

•  The  problem  is  restricted  to  one  spatial  dimension. 

•  Charge  neutrality  is  maintained. 
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•  The  neutral  density  and  velocity  profiles  are  prescribed  by  a  step  function.  In 
front  of  the  shock  the  fluid  is  static,  while  behind  the  shock  the  neutral  velocity 
is  given  by  u. 

The  complete,  nonmagnetized  ion  and  electron  momentum  equations  are  given 
below,  where  V  is  the  fluid  velocity,  v  is  the  collision  frequency,  and  the  subscripts, 
n,  i,  and  e  represent  neutrals,  ions,  and  electrons,  respectively. 


^(PiVi)  +  -  Ve )  -  Pi”in{yi  -  yn)  (24) 


~^(Peye)  +  K^(Pe^e)  =  ~ V-Pe  ~  “ - Pe^ei(ye  ~  PeUen( K  K.)  (25) 

After  applying  the  above  constraints,  the  electron  momentum  equation  reduces  to 


kB 

me 

and  the  electric  field  can 

E 


^  ^  -  PeUen{yi  ~Vn)  =  0 

OX  ““ 


me 


be  approximated  as 


fc b  9(peTe)  —  meven  —  V  ) 
pee  dx  e 


(26) 


(27) 


The  two-fluid  approximation  is  completed  after  substituting  the  electric  field  into 
Equation  24,  and  assuming  quasi-neutrality,  in  which  ne~  up. 


+  V, ^(PiVi)  -  ~ 


ksdip^l  _  +  j)  _  Vn)  (28) 

m  ox  m  uin 
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The  collision  frequency  is  given  by  Gombosi  (17:99)  as 


V 12  =  2(7 12712 


(29) 


where  the  subscript  1  represents  the  particle  that  collides,  subscript  2  represents 
the  particle  that  is  impacted,  and  a  is  the  collision  cross  section.  When  applying 
Equation  29  to  ion  momentum  equation,  the  mass  and  collision  frequency  ratio  term 
is  much  less  than  unity  and,  therefore,  neglected.  Next,  the  one-dimensional  ion 
pressure  gradient  is  isothermally  expanded  using  the  ideal  gas  law  (Equation  4). 
Then,  the  ion  temperature  is  dropped  since  Te  »  T,,  the  convection  derivative  is 
expanded,  and  the  non-steady  term  is  dropped.  Finally,  recall  the  ion-acoustic 
velocity,  Equation  23,  and  the  charged  component  momentum  becomes: 

+  V^)  =  -(vi  +  V?)||  -  wUVi  -  vn)  (30) 


The  hydrodynamic  fluid  velocity,  V,  is  never  greater  than  c,  the  shock  velocity; 
therefore  V  <  Ma.  The  ratio  V?/v2ia  is  approximately  7 M2Tj/Te.  Since  electron 
temperatures  are  usually  two  orders  of  magnitude  greater  than  ion  temperatures, 
and  if  the  Mach  number  remains  near  unity,  the  ratio  is  much  less  than  unity  and 
the  ion  velocity  is  neglected.  Therefore,  Equation  30  becomes 


m  w 

dt  1  dx 


Di  dx 


ViniVi  Fn) 


(31) 


Since  Avramenko  uses  a  prescribed  solution  for  the  neutral  velocity  and  if 
the  upstream  ion-neutral  collision  frequency,  vin,  can  be  treated  as  a  parameter, 
then  there  are  two  unknowns,  Vi  and  p{.  The  continuity  equations  for  the  neutral 
and  charged  components  round  out  the  required  number  of  equations  to  solve  the 
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problem: 


f+!<^=°  (32) 

^  +  |(,„K)  =  °  (33) 

Even  though  Equation  33  is  a  third  equation,  it  will  be  used  to  establish  boundary 
conditions. 

The  equations  are  solved  in  the  shock-fixed  frame;  therefore,  it  is  necessary  to 
introduce  the  transformation  equations.  The  independent  variable  in  this  frame  is 
£  =  x  —  ct,  where  c  is  the  shock  velocity  in  the  “laboratory  frame.  In  the  shock 
frame,  the  shock  remains  fixed  at  £  =  0,  while  the  upstream  fluid  is  moving  at 
V  =  —  c.  The  differential  transformations  are: 

8_  _  d_d£  d_ 

9t~  d£dt~  °d 

The  velocity  transformations  are: 

Vi  yi  =  c  -  Vi, 

where  yn  and  yi  are  the  neutral  and  ion  velocities  in  the  shock  frame,  respectively. 

Now  both  ion  and  neutral  continuity  equations  and  the  momentum  equation 
can  be  transformed.  Since  £  incorporates  both  variables  x  and  t,  the  partial  deriva¬ 
tives  are  also  transformed  into  total  derivatives. 

-cf +|(ft(c'!,i))  =  0 

£L  =  Vi!  = - —  (36) 

Pio  yi  c- Vi 


d_  _  d_d^_  _  d_ 
dx  d£  dx  d( 


(34) 


Vn  ^  yn  —  C  Vn 


(35) 
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Here,  pi0  is  the  undisturbed  upstream  ion  density  and  the  ion  velocity  far  upstream 
of  the  shock,  Vi0,  is  assumed  to  be  zero.  The  neutral  continuity  equation  transforms 
similarly,  with  the  exception  that  neutral  velocity  and  density  profiles  are  prescribed 
by  step  functions  at  the  shock. 

-C^  +  -^(Pn(C“2/n))  =  0 

d£  d£ 

Pn  1  _  Vn 0  _  C  (37) 

PnO  Vnl  C  ~  U 

Here,  pni/pn0  is  the  neutral  density  ratio  across  the  shock  and  u  is  constant  neutral 
velocity  behind  the  shock.  Equation  31  becomes 

d£  Pi 

«d)  +  ±d-B.-vUyi-yn)  =  0  (38) 

d£  ^  Pi 

The  shock-frame  portion  of  Equation  36  is  now  inserted  into  the  ion  momentum 
equation: 

Y  -  V2ia  ln(“))  -  ViniVi  -Vn)=  0  (39) 

Due  to  the  neutral  step  functions,  this  differential  equation  must  be  solved  in 
two  parts:  for  the  post-shock  region,  £  <  0,  and  the  upstream  region,  £  >  0,  with 
continuity  of  both  V*  and  p,  maintained  at  £  =  0.  The  post-shock  region  is  solved 
first.  The  key  to  the  solution  is  the  assumption  that  there  are  no  forces  acting  on 
the  fluid  behind  the  shock;  therefore,  d^/d£  =  0.  In  the  laboratory  frame,  the 
neutral  velocity  is  a  constant  value  of  u  in  this  region;  therefore,  yr  =  c  —  Vi  and 
yn  =  c-u.  After  introducing  these  relations,  expressing  the  ratio  Vi/c  as  V ,  and 
simplifying  under  the  assumption  that  V, 2  Ma  <  1,  Equation  39  can  be  rearranged 
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as  such: 


dV_  _  cum  (H  _  y-)(i  _  fr)  =  Q  (40) 

d£v2iac 

The  two  possible  solutions  are  V  =  1  or  V  =  u/c.  The  first  solution  means  that 
Vi  =  c ;  this  is  only  possible  for  shocks  moving  into  a  vacuum,  where  pn0  =  0.  Since 
a  general  solution  is  sought,  then  V  =  u/c  and  Vi  =  u.  Applying  the  neutral 
continuity  equation,  Equation  37,  yields  the  following  solution: 


V  =  u  =  c(l  -  ^0)  (41) 

Pn  1 

for  the  post-shock  ion  velocity  in  the  laboratory  frame.  Equation  36  provides  the 
post-shock  ion  density: 


c  c  _  Pni 

PiO~  77  ~~  Pi0  n  0,  ~  ^i0  n 
C~  Vi  c  —  u  pn0 


(42) 


It  is  now  possible  to  obtain  the  solution  upstream  from  the  shock.  In  the 
laboratory  frame,  the  neutral  velocity  is  assumed  to  be  zero  in  this  region;  therefore, 
Ui  =  c  — Vi  and  yn  =  c.  After  substituting  these,  V  =  V/c,  and  applying  the 
assumption  V?/v?a  1,  Equation  39  becomes 


1  dV  _  cuin 
V(1  -  V)  Mi  ~  v\a 


(43) 


a 

Integrating  from  £  =  0  where,  from  continuity  with  the  post-shock  region,  V  =  u/c, 
to  some  positive  £  and  V,  yields  the  following, 


Ju/c  V(1  -  V)  -  Jo 


V) 

V 


vt 


l-d£ 


!  +  (£-!)  exp(^£)  1  +  (tsizt)  exp(^£) 


(44) 
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Avramenko  introduced  the  following  quantities  to  simplify  the  solution: 


n  =  H—  - 1).  U  =  -%  <45) 

Pn 0  CU\n' 

where  f0  is  the  characteristic  width  of  the  precursor  and  v$  is  taken  to  be  ambient 
upstream  ion-neutral  collision  frequency,  given  by  Equation  29.  The  solution  of  the 
upstream  ion  velocity  in  the  laboratory  frame  is 

V-  = _ _ _  (46) 

1  +  exp(£  -  n) 

Calling  upon  Equation  36,  the  upstream  ion  density  is  given  by: 

Pi  =  Pio—^TF  =  Ao(l  +  “  |-)) 

c  —  Vi  So 


Now,  examine  the  solution  for  a  shock  where  M  —  2  for  argon,  where  7  5/3 

and  a  =  10-6.  Equation  17  yields  the  following  neutral  ratios:  p2/pi  =  2.285  and 
u2/ui  =  0.438.  The  ambient  upstream  conditions1  of  Ti<n  =  300  K,  Te  =  2  eV (= 
23200  K) ,  Pn 0  =  30  torr,  and  ain  =  4  x  10-21  m2,  yield  a  sound  speed  of  c  =  323  m  /  s, 
a  shock  speed  of  c  =  645  m  /  s,  and  a  precursor  width  of  £0  =  3.4  mm.  The  solution 
of  the  charged  component  velocity  is  shown  in  Figure  6,  and  the  density  is  shown  in 
Figure  7. 

Recall  Equation  27,  the  electric  field  approximation.  Hilbun  has  shown  that 
the  collision  term  is  negligible  (22:75)  and  the  field  approximation  is  easily  converted 
to: 


fcg  djfijTe) 
epi 


(48) 


1  Unless  otherwise  indicated,  these  values  will  serve  as  the  ambient  upstream  conditions  for  the 
remainder  of  the  present  research. 
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Figure  6  Steady  state  analytical  solution  of  charged  component  velocity  for  a  Mach 
2  shock  in  argon.  Velocity  is  normalized  to  the  shock  velocity,  c.  Dis¬ 
tance  is  normalized  by  £0  =  3.4  mm. 
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Figure  7  Steady  state  analytical  solution  of  charged  component  density  for  a  Mach 
2  shock  in  argon.  Density  is  normalized  to  the  upstream  charged  com¬ 
ponent  density.  Distance  is  normalized  by  £0  =  3.4mm. 
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Figure  8  Steady  state  analytical  approximation  of  the  electric  field  in  the  vicinity 
of  a  Mach  2  shock  in  argon.  Distance  is  normalized  by  f 0  =  3.4  mm. 

In  the  post-shock  region,  the  field  approximation  is  zero  since  p{  is  constant.  In  the 
fax  upstream  regions  the  field  is  zero  because  p±  =  pi0.  Just  ahead  of  the  shock,  the 
electric  field  becomes 


ksTe  ^  _  PiO  ^  Vj 

eCo  Pi  <0  C 


(49) 


Recalling  Equation  36,  the  electric  field  approximation  has  the  same  functional  form 
as  the  charged  component  velocity  in  the  upstream  region,  as  shown  in  Figure  8. 
Although,  in  reality,  since  the  electric  field  drives  the  charged  particles,  the  charged 
component  velocity  follows  the  form  of  the  electric  field.  The  potential  drop,  as 
determined  by  Equation  3,  experienced  by  a  charged  particle  in  this  electric  field  is 
about  —1.5  V. 


3.1.2  Steady  State  Numerical  Analysis  .  Hilbun  (22:69)  initially  obtained 
estimates  of  the  flow  field  variables  and  electric  field  by  an  extension  of  Avramenko’s 
steady-state  approach.  He  refined  the  neutral  velocity  and  density  profiles  at  the 
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shock  front  with  a  more  physical  profile-a  Fermi  function,  given  as  follows 


(50) 

(51) 


where  £  defines  the  shock  width,  which  is  usually  on  the  order  of  14  neutral-neutral 
mean 'free  paths.  This  restriction,  however,  still  prevents  the  charged  compo¬ 
nent  from  affecting  the  neutral  component.  He  also  removed  the  assumption  that 
V?/vfa  <  1,  which  allows  for  stronger  shocks.  With  these  assumptions,  Equation 
39  becomes 


dyi  _  Vinjyi  -  yn{Q)yj 
_  Vi  -  Via 


This  equation  yields  smooth  solutions  of  the  charged  component  velocity  and  density. 
From  these,  Equation  48  yields  a  smooth  electric  field  estimate  through  the  shock. 
The  continuous  electric  field  allows  the  derivation  of  the  net  charge  density  in  the 
vicinity  of  the  neutral  shock  front  by  way  of  Poisson’s  equation,  Equation  2. 

These  extensions  require  a  numerical  solution.  Equation  52  was  solved  using  a 
fourth-order  Runge-Kutta  routine.  The  collision  frequency  is  given  by  Equation  29. 
Figure  9  depicts  the  charged  component  velocity.  The  charged  component  density 
is  shown  in  Figure  10.  Note  that  the  numerical  and  analytical  solutions  show  only 
minor  deviations  for  the  parameters  used.  The  electric  field  is  depicted  in  Figure  11. 
A  charged  particle  passing  through  this  field  would  experience  a  potential  drop  of 
—1.5  V,  which  is  in  good  agreement  with  the  analytical  treatment.  Finally,  Figure 
12  illustrates  the  separation  of  charge  that  occurs  at  the  shock  front.  The  electrons 
are  broadly  distributed  through  the  precursor  due  to  their  high  mobility,  while  the 
ions  are  compacted  at  the  shock  front  due  to  the  strong  electric  field  and  their  low 
mobility. 
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Figure  9  Steady  state  numerical  solution  of  charged  component  velocity  (broad 
dashes)  for  a  Mach  2  shock  in  argon.  Analytic  solution  (solid)  and  neu¬ 
tral  velocity  (short  dashes)  shown  for  comparison.  Velocity  is  normalized 
to  the  shock  velocity,  c.  Distance  is  normalized  by  £0  =  3.4  mm. 
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Figure  10  Steady  state  numerical  solution  of  charged  component  density  (broad 
dashes)  for  a  Mach  2  shock  in  argon.  Analytic  solution  (solid)  and 
neutral  density  (short  dashes)  shown  for  comparison.  Densities  are  nor¬ 
malized  to  upstream  values.  Distance  is  normalized  by  £0  =  3.4  mm. 
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Figure  11  Steady  state  numerical  solution  of  the  electric  field  in  the  vicinity  of  a 
Mach  2  shock  in  argon.  Distance  is  normalized  by  £0  =  3.4  mm. 
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Figure  12  Steady  state  numerical  approximation  of  the  net  charge  density  in  the 
vicinity  of  a  Mach  2  shock  in  argon.  Distance  is  normalized  by  £0  = 
3.4  mm.  Electronic  double  layer  is  clearly  visible. 
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3.2  Time- dependent 

Hilbun’s  development  continues  with  a  time-dependent  approach.  In  his  treat¬ 
ment  of  shock  propagation,  he  examined  three  facets  of  weakly  ionized  gases-a  vi- 
brationally  excited  gas,  a  hot  gas,  and  a  plasma-and  developed  numerical  solutions 
for  each  case  (22:25).  In  the  case  of  the  vibrationally  excited  gas,  Hilbun  assumed 
the  gas  is  a  non-equilibrium  store  of  vibrational  energy.  The  objective  was  to  de¬ 
termine  the  amount  of  vibrational  energy  that  is  transferred  to  translational  energy 
as  the  shock  passes  through  it.  For  the  hot  gas  approximation,  Hilbun  used  a 
two  dimensional  fluid  dynamics  code  to  investigate  the  modifications  on  the  shock 
propagation  parameters  as  it  passes  through  thermal  inhomogeneities.  The  present 
research  focuses  on  the  plasma  aspect  of  the  weakly  ionized  gas. 

3.2.1  One  dimensional,  two- fluid  approximation.  Hilbun’s  approach  to 
modelling  the  medium  is  also  a  two-fluid  approximation-neutrals  and  charged  par¬ 
ticles,  where  the  electrons  are  assumed  to  be  steady-state.  Therefore,  the  electric 
field  can  be  derived  from  the  steady  state  electron  momentum  equation,  and  is  given 
by  Equation  48.  The  treatments  reviewed  in  the  previous  section  assume  static 
neutral  density,  velocity,  and  temperature  profiles;  however,  the  present  approach 
allows  for  modification  of  the  neutral  flow.  Interaction  between  the  fluids  occurs 
through  energy  and  momentum  coupling  and  there  are  no  production  or  loss  mech¬ 
anisms.  Viscosity  and  thermal  conductivity  are  assumed  to  be  negligible,  which 
allows  Euler’s  equations  to  be  used  instead  of  the  complex  Navier-Stokes  equations. 
Therefore,  the  neutral  profiles,  and  those  of  the  charged  component,  are  determined 
by  solving  Euler’s  equations. 

Hilbun’s  time-dependent  equations  take  on  the  following  form  (22:79): 


dU  dF_ 
dt  dx 


=  Si  +  s2  +  sz 


(53) 
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where  the  terms  on  the  left  are  the  derivatives  of  the  conserved  variables-mass, 
momentum,  and  energy-and  the  terms  on  the  right  are  called  source  terms-sources 
of  mass,  momentum,  and  energy.  The  specific  vectors  are  defined  below.  The 
vector  of  conserved  variables,  U,  is  defined  as 


The  source  terms  for  the  Euler  equations  and  species  coupling  are  given  by  the 
following: 
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where  Pjk  represents  the  momentum  gained  by  species  j  at  the  expense  of  species  k, 
and  Qjk  represents  a  similar  transfer  of  energy  between  the  two  species.  The  mo¬ 
mentum  transfer  between  ions  and  electrons  is  zero  by  definition  since  their  velocities 
are  assumed  to  be  equal  in  the  two-fluid  approximation.  Jaffrin  (23:611)  gives  the 
specific  momentum  and  energy  coupling  terms  as 

Pni  =  7finni(Jin{Vi  -  Vn)^kB(Tn  +  T,)  (57) 

Pne  =  |nnne(7ew(Vri  -  Vn)^2me^B—  (58) 


Qni  =  2nnrii(Tin \/——  kB (Tn  4-  Ti){kB(Tn  -  Ti )  +  -m(Vi  — 14) (V*  + 14)}  (59) 

V  7rm  o 


Qm  =  Snnne*„]fP^{^(Te  -  T„)  +  -  V.)(V.  +  ^V,)}  (60) 


Qei  =  8riineaei—(Ti  -  Te ) 
m 


‘ZmekiiTe 

7T 


The  electric  field  source  term  is  given  by 


S3  = 


kBTe  d£i 

m  dx 


0 

0 

0 

0 

1 


(61) 


(62) 
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where  the  electric  field,  E,  is  approximated  by  Equation  48. 

Hilbun’s  numerical  techniques  involve  a  finite  difference,  second  order  accu¬ 
rate  MacCormick  algorithm  with  flux  corrected  transport  (FCT).  The  MacCormick 
method  is  widely  used  in  solving  fluid  flow  equations  (3:102).  It  is  a  predictor- 
corrector  method,  in  which  the  predictor  term  uses  forward  differencing  and  the 
corrector  term  uses  backwards  differencing.  This  differencing  can  be  reversed  for 
cases  involving  moving  discontinuities,  such  as  shocks.  The  FCT  algorithm  rids 
the  solution  of  oscillations  resulting  from  second  order  algorithms.  The  algorithm 
advances  the  solution  by  a  time  step  in  a  two-step  process  (38:84).  In  each  half-step, 
a  diffusive  flux  is  applied  to  the  solution  to  introduce  numerical  diffusion  to  ensure 
stability  and  monotonicity,  and  then  an  anti-diffusive  flux  is  applied  to  eliminate 
excessive  numerical  flux. 

In  order  to  compare  the  solution  obtained  with  Hilbun’s  plasma  code  to  Avra¬ 
menko’s  solution,  a  Mach  2  test  case  was  completed.  The  initial  conditions  Pj Pi  = 
22.2  and  T4/T1  =  2.0  generate  a  Mach  2  shock  in  a  shock  tube,  as  determined  by 
Equation  19.  Hilbun’s  time-dependent  code  accurately  solves  the  Riemann  problem, 
as  illustrated  in  Figure  13  (compare  to  Figure  4),  though  there  were  no  modifications 
to  neutral  shock  structure  and  propagation.  The  time-dependent  solution  matches 
Avramenko’s  results  in  that  a  charged  precursor  appears  before  the  shock,  as  shown 
in  Figure  14.  The  solution  also  yields  an  electric  field  near  the  shock,  as  shown  in 
Figure  15.  Though  the  field  may  be  broader  and  less  intense  than  the  field  obtained 
with  the  steady-state  approach,  a  charged  particle  passing  through  this  field  would 
also  experience  a  potential  drop  of  — 1.5  V-  Net  charge  density  is  approximated  by 
Equation  2  and  is  illustrated  in  Figure  16;  the  EDL  is  clearly  visible,  though  the  ionic 
portion  of  the  EDL  is  far  less  concentrated  than  in  the  steady  state  case  (compare  to 
Figure  12).  Since  this  approach  assumes  charge  neutrality,  the  globally  integrated 
charge  density  should  sum  to  zero.  Letting  the  numerical  error  be  measured  with 
the  respect  to  the  maximum  charge  density  in  Figure  16,  charge  neutrality  is  mam- 
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Figure  13  Hilbun’s  one-dimensional  plasma  code  solution  of  neutral  pressure  for  a 
Mach  2  shock  in  argon.  Shock  and  expansion  waves  are  clearly  visible. 

tained  to  within  98.8  percent.  The  magnitude  of  this  error  is  due  in  part  to  the 
coarseness  of  the  spatial  grid  used  in  the  numerical  solution. 

There  are  two  minor  concerns  when  using  Hilbun’s  plasma  code.  First,  mo¬ 
mentum  and  energy  coupling  between  the  heavy  particles  and  electrons  must  be 
neglected.  Coupling  is  afforded  in  the  code  via  Equations  58,  60,  and  61.  Since  the 
electron  temperature  is  assumed  constant,  the  electron  coupling  terms  provide  an 
unlimited  source  of  energy  for  the  heavy  particles  as  Tn  and  T*  equilibrate  with  Te. 
Therefore,  in  two-fluid  calculations,  it  is  best  to  neglect  electron  coupling  by  setting 
Oei  and  (Ten  to  zero.  Secondly,  using  an  average,  order  of  magnitude  approximation 
of  the  ion-neutral  collision  cross  section,  ain  =  10~19  m2,  causes  numerical  instabil¬ 
ities.  The  code  requires  a  cross  section  that  is  25  times  smaller  than  this  value 
(22:89).  This  is  most  likely  due  to  the  fact  that  the  average  time  step  in  the  calcula¬ 
tion,  At  Ax/ Via ,  is  larger  than  the  ion-neutral  collision  time,  Tjn  ~  1  /■v/2nn<7tT,Vj. 
Assuming  an  average  o jn  of  10-19  m2  and  test  conditions  of  Te  =  2  eV,  Tt  =  300  K, 
and  Pn0  =  30  torr,  the  ratio  At/rin  is  on  the  order  of  102.  Nominal  values  of  ain 
could  be  used  if  the  computational  time  step  is  reduced  by  two  orders  of  magnitude. 
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Figure  14  Solution  of  charged  component  density  (dashed)  for  a  Mach  2  shock  in 
argon.  Neutral  density  (solid)  shown  for  comparison.  Densities  are 
normalized  by  upstream  values. 


Figure  15  Approximation  of  electric  field  for  a  Mach  2  shock  in  argon. 
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Figure  16  Approximation  of  net  charge  density  for  a  Mach  2  shock  in  argon. 

Unfortunately,  the  time  and  computational  restraints  levied  upon  this  research  effort 
precluded  the  use  of  these  more  accurate  calculations. 
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IV.  Lifted  Restrictions 

According  to  Hilbun  (22:96),  in  order  for  ion-acoustic  wave  damping  to  affect  the 
neutral  flow,  two  conditions  must  be  met:  1)  The  ion-acoustic  wave  energy  density 
must  be  an  appreciable  fraction  of  the  neutral  energy  density,  and  2)  the  time  scale 
for  the  energy  transfer  from  the  ion- acoustic  wave  to  the  neutral  flow  must  be  less 
than  the  transit  time  for  a  neutral  particle  to  pass  through  the  charged  component 
precursor.  Hilbun  derives  the  threshold  fractional  ionization,  aia,  that  meets  these 
two  criteria:  ctja  ^  5x  10— 3.  Figure  17  illustrates  the  effect  of  ion-neutral  momentum 
coupling.  Various  fractional  ionizations  were  tested  using  the  same  initial  conditions, 
PA/P1  =  10.0  and  T4/T1  =  1.25  1.  The  horizontal  line  represents  shock  speed  that 
is  unaffected  by  the  presence  of  the  plasma,  and  the  vertical  line  represents  Hilbun  s 
value  of  aia.  Shock  acceleration  due  to  ion-acoustic  damping  is  clearly  visible  for 
Hilbun’s  value  of  aia,  but  the  results  of  his  time-dependent  solution  seem  to  favor 
ocia  ~  1  x  10-3  as  the  threshold. 

According  to  Adamovich,  et  al.  (1:816),  and  Jaffrin  (23:610),  gaseous  media 
with  a  ~  10-3  are  generally  classified  as  partially  ionized  gases,  which  are  not 
considered  within  the  scope  of  this  research.  However,  if  the  value  of  a  of  a  weakly 
ionized  gas  were  to  increase  to  partially  ionized  values  in  the  vicinity  of  the  shock, 
then  ion-acoustic  wave  damping  may  be  effective  in  altering  shock  structure  and 
propagation.  A  possible  mechanism  for  varying  the  fractional  ionization  at  the  shock 
front,  as,  is  additional  ionization.  Recall  the  modified  ion  momemtum  equation  from 
the  previous  chapter  (Equation  24): 

f  WO  +  V'~uvi)  =  -vp,  +  ^E-  »(V,  -  vi) 

1These  initial  conditions  will  be  used  for  the  remainder  of  the  work.  Unless  otherwise  stated, 
all  plasma  code  solutions  assumed  a  shock  tube  length  of  lm,  with  the  diaphragm  located  at 
x  —  0.4  m,  and  all  shocks  propagated  to  x  —  0.9  m. 
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Figure  17  Neutral  shock  propagation  speed  vs.  fractional  ionization,  a.  The  ver¬ 
tical  line  indicates  the  threshold  fractional  ionization,  aia,  for  effective 
ion-acoustic  wave  damping  as  derived  by  Hilbun.  The  horizontal  line 
indicates  nominal  shock  velocity  in  an  unionized  gas. 

Both  the  pressure  gradient  and  electric  field  forces  act  in  the  direction  of  shock 
propagation.  It  was  made  plain  by  Figure  15  that  there  is  a  positive  electric  field 
at  the  shock  front.  These  forces  work  to  increase  the  ion  velocity  above  the  neutral 
velocity.  As  the  ion  velocity  increases,  however,  collisional  coupling  with  the  neutrals 
restricts  the  forward  movement  of  the  ions.  Because  of  this,  a  large  portion  of  the 
ions  generated  at  the  shock  front  will  tend  to  remain  there.  As  the  ion  density 
builds  at  the  shock,  it  could  reach  a  level  conducive  to  momentum  transfer  to  the 
neutrals  via  ion-acoustic  wave  damping. 

In  Hilbun’s  conclusion  on  plasma  effects  (22:96),  he  states  that  in  his  time- 
dependent  solution,  “Under  plasma  conditions  typically  encountered  in  glow  dis¬ 
charges,  the  charged  particles  are  observed  to  have  practically  no  influence  on  the 
density,  velocity,  and  temperature  of  the  neutral  component  in  the  parameter  space 
investigated.  The  neutral  shock  velocity  also  remained  unaffected  by  the  plasma 
component  under  these  conditions.”  The  field  strength  was  comparable  to  the  axial 
field  of  a  glow  discharge,  but  was  of  such  small  spatial  extent  that  electrons  were 
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not  significantly  affected  (22:75).  Hilbun’s  argument  against  additional  ionization  is 
this:  The  peak  electric  field  at  the  shock  front  is  of  the  same  order  of  magnitude  as 
the  axial  electric  field  typically  found  in  a  glow  discharge  in  argon  at  30torr,  about 
50  V /cm.  However,  the  field  is  localized  to  a  small  region  around  the  shock  front; 
therefore,  the  total  potential  drop  is  only  on  a  few  Volts.  With  electron  thermal 
energies  of  about  2  eV,  this  potential  drop  falls  short  of  ionization  potential  of  argon, 
IP  =  15.8  eV  (40:389),  and  is  unable  to  significantly  increase  ionization  at  the  shock 

front. 

The  logical  extrapolation  of  this  argument  is  that  stronger  electric  fields  and 
greater  precursor  widths  should  lead  to  greater  potential  drops  for  electrons.  These 
conditions  could  be  generated  by  elevated  electron  temperatures  at  the  shock  front. 
The  resulting  potential  drops  could  approach  the  ionization  potential  of  argon. 
Adopting  a  quasi-kinetic  approach,  for  the  given  average  thermal  energy  of  the  elec¬ 
trons  (about  2eV)  a  small  fraction  of  the  electrons  will  have  energies  greater  than 
IP.  This  fraction,  albeit  small,  is  3.7  x  10~4  and  is  given  by  the  Boltzmann  factor, 
exp(-IP/kBTe).  If  the  ion  density  builds  over  time,  ion-neutral  momentum  cou¬ 
pling  could  begin  to  play  a  significant  role  in  shock  structure  and  propagation  via 
ion-acoustic  wave  damping,  which  was  mentioned  in  Section  2.2.2. 

It  may  seem  intuitous  that  larger  amounts  of  plasma  would  generate  larger 
electric  fields  due  to  the  greater  numbers  of  charged  particles  that  constitute  the 
charge  separation.  This,  however,  is  not  the  case.  From  the  previous  chapter,  the 
electric  field  approximation  is  given  by 

kB  d(PjTe)  _  _/cb,7; dn i  +  (63) 

The  electric  field  is  directly  related  to  the  charged  component  pressure  gradient,  in 
which  the  electron  thermal  energy  and  its  spatial  gradient,  and  the  spatial  gradient 
of  the  ion  density  are  all  factors.  More  energetic  electrons  and  stronger  shocks 
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Figure  18  Comparison  of  electric  fields  for  a.  =  10  2  (solid)  and  a.  —  10  6  (dashed). 

Charged  particles  in  each  of  these  fields  experience  a  potential  drop  of 
only  -IV. 

will  generate  stronger  electric  fields.  The  field  is  independent  of  the  actual  ion 
density,  however.  As  a  increases,  the  electric  field  and  potential  remain  relatively 
unchanged,  which  is  the  case  in  Figure  18.  For  both  values  of  ol  represented,  the 
associated  electric  fields  both  produced  potential  drops  on  the  order  of  —1 V*  Also, 
as  a.  increases,  the  charge  separation  and  net  charge  density  are  relatively  unchanged, 
as  seen  in  Figure  19.  Yet,  in  Figure  17,  the  cases  in  which  a>lx  10-3,  the  shocks 
propagated  at  velocities  in  excess  of  the  unionized  Riemann  case.  This  indicates  that 
ion-neutral  momentum  coupling  is  the  critical  factor  in  the  modification  of  neutral 
shock  structure  and  propagation.  The  required  increase  in  ion  density,  however, 
can  not  come  about  unless  the  electric  potential  is  of  sufficient  strength  to  enhance 
electron  impact  ionization.  The  field  strength,  of  course,  is  a  function  of  the  electron 
energy  at  the  shock  front. 

It  is  the  intent  of  this  research  to  extend  the  investigation  of  the  structure  and 
propagation  of  shocks  in  weakly  ionized  gas  as  the  EDL  is  freely  affected  by  variable 
electron  temperature  and  additional  electron  impact  ionization  at  the  shock.  Greater 
ion-neutral  coupling,  assisted  by  a  stronger  electric  field  and  potential,  could  provide 
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Figure  19  Comparison  of  net  charge  density  profiles  for  a  =  10-2  (solid)  and 
a  —  10-6  (dashed). 

the  mpa.ns  for  shock  modification  in  weakly  ionized  gases.  A  spatially-dependent 
electron  temperature  will  enhance  the  strength  and  breadth  of  the  electric  field,  as 
evident  in  Equation  63.  The  enhancements  in  electron  temperature  will  sustain 
ionization  at  the  shock  front  by  increasing  the  energies  of  free  electrons.  This 
chapter  investigates  these  mechanisms.  Each  of  the  following  sections  are  divided 
into  the  development  of  the  lifted  restriction,  a  description  of  the  results,  and  an 
analysis  of  the  results. 

4-1  Variable  electron  temperature 

4-1.1  Development.  Recent  research  by  Sirghi,  et  al.,  (36),  indicates  that 
electron  temperatures  increase  in  the  vicinity  of  a  strong  electric  field.  The  goal 
of  their  research  was  to  determine  the  distribution  function  of  electrons  at  sharp 
changes  in  diameter  in  glow  discharge  tubes.  In  their  working  media  of  helium 
at  1  torr,  large  electric  fields  were  generated  at  these  constrictions,  with  values  two 
to  four  times  the  ambient  homogeneous  positive  column  electric  field.  With  an 
homogeneous  field  in  DC  positive  column  of  about  4  V  /  cm,  the  induced  field  expe¬ 
rienced  a  maximum  of  approximately  11 V  /  cm  ahead  of  the  constriction,  and  had  a 
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Figure  20  Electric  field  vs.  distance.  Rise  in  field  was  generated  by  a  constriction 
in  a  glow  discharge  tube.  Field  structure  is  similar  to  that  experienced 
at  a  shock  in  weakly  ionized  gas  (Ref:  (36).  Note  that  geometry  has 
been  reversed). 

width  of  approximately  1  cm.  This  field  structure  is  illustrated  in  Figure  20,  where 
the  mouth  of  the  constriction  resides  at  x  =  0,  and  the  constricted  region  is  where 
x  <  0.  This  induced  field  is  not  unlike  the  electric  field  generated  by  the  EDL  at  a 
shock  front  in  a  weakly  ionized  gas.  In  the  region  of  elevated  electric  field,  ambient 
electron  energies  of  approximately  5  eV  were  nearly  doubled  as  seen  in  Figure  21. 
Sirghi  attributed  the  increased  electron  energy  to  the  induced  electric  field.  It  is 
clear  from  Equations  58,  60,  and  61  that  momentum  and  energy  coupling  with  the 
heavy  particles  cannot  account  for  the  distributions,  since  this  would  only  heat  the 
heavy  particles  and  lower  the  energies  of  the  electrons  due  to  equilibration.  Even 
though  the  research  cited  here  involves  weakly  ionized  helium,  similar  phenomena 
should  be  expected  in  argon  as  well. 

The  restriction  of  constant  electron  temperature  was  lifted  by  incorporating 
a  spatially-dependent  region  of  electron  heating  in  the  vicinity  of  the  neutral  shock 
front.  This  region  was  similar  in  form  to  the  electron  thermal  layer  evident  in  Figure 
21  and,  for  simplicity,  took  the  form  of  a  shock-centered  Gaussian  rise  in  electron 
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Figure  21  Electron  energy  vs.  distance.  Rise  in  energy  was  generated  by  a  con¬ 
striction  in  a  glow  discharge  tube  (Ref:  (36).  Note  that  geometry  has 
been  reversed). 

temperature  over  and  above  the  baseline  electron  temperature,  Te q: 


Te{0  =Te o  +  ATeexp(-(-^)2) 


(64) 


where  ATe  is  the  increase  in  electron  temperature,  usually  a  fraction  of  Te o,  and  6 
is  the  number  of  upstream  mean  free  paths,  A,  that  specify  the  approximate  half¬ 
maximum  width  of  the  Gaussian  curve.  A Te  and  6  serve  as  the  coordinates  of  the 
parameter  space  investigated.  The  mean  free  path  is  given  by  (17:95)  as: 


A  = 


1 


VZnrtoCTin 


(65) 


where  nn o  is  the  upstream  neutral  number  density  and  a in  is  the  ion-neutral  collision 
cross  section. 

The  central  assumption  of  the  two-fluid  approximation  was  that  of  steady- 
state  electron  momentum.  It  might  seem  that  the  incorporation  of  a  spatially- 
dependent  electron  temperature  would  violate  this  assumption.  It  is  important  to 
note  that  the  assumption  is  not  violated.  Consider  Equation  53  and  its  constituent 
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vectors  as  a  basis.  On  the  source  side  of  the  equation,  if  the  collisional  coupling 
of  the  electrons  is  ignored,  then  the  momentum  is  maintained  by  a  balance  between 
the  pressure  gradient  and  the  force  imparted  by  the  electric  field.  The  variable 
temperature  introduces  a  temperature  gradient  term  in  the  pressure  gradient,  which 
is  compensated  also  by  a  temperature  gradient  term  in  the  electric  field.  On  the 
conserved  variable  side,  changes  in  density  must  be  maintained  by  changes  in  fluid 
velocity.  Furthermore,  the  conservation  of  energy  is  maintained  as  well.  For  the 
source  side  of  Equation  53,  a  balance  between  the  pressure  gradient  and  the  electric 
field  term  similar  to  that  described  above  holds  here  as  well.  On  the  left  side, 
a  balance  between  electron  temperature,  density,  and  fluid  velocity  is  maintained. 
Of  note,  however,  is  the  fact  that  the  strengths  and  spatial  extents  of  the  electron 
thermal  layer  are  based  on  empirical  results  rather  than  a  solution  of  the  electron 
continuity,  momentum,  and  energy  equations.  Temperature  increases  and  gradients 
that  are  too  large  may  lead  to  nonphysical  solutions  of  density  and  velocity. 

The  spatially-dependent  electron  temperature  was  easily  incorporated  into 
Hilbun’s  two-fluid  plasma  code.  A  simple  subroutine  employing  Equation  64  was 
implemented.  The  value  of  A  was  calculated  in  the  code  in  accordance  with  Equation 
65,  the  parameters  A Te  and  <5  were  user-supplied,  and  the  shock  location,  £  =  0,  was 
determined  in  the  code  for  each  time  step.  In  addition  to  the  Te(£)  modification, 
the  electric  field  source  term,  S3,  was  modified  to  account  for  the  variable  electron 
temperature: 


Ss  = 

m  ox 


0 

0 

0 

0 

1 


(66) 
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Figure  22  Normalized  ion  density  profiles  resulting  from  variable  electron  temper¬ 
ature.  Constant  Te  =  2 eV  (solid  line),  A Te  =  1  eV  (short  dashes),  and 
A Te  =  2eV  (broad  dashes). 

4.1.2  Results.  Incorporating  the  spatial  variation  of  the  electron  tempera¬ 
ture  led  to  three  effects  on  the  charged  component  flow  parameters  and  electric  field 
structure.  The  first  was  an  increase  in  the  charged  component  precursor  width,  £0. 
This  increase  in  precursor  width  is  illustrated  in  Figure  22,  where  the  solid  line  rep¬ 
resents  the  normalized  ion  density  of  the  constant  Te  =  2  eV  case,  the  short  dashes 
represent  A Te  =  1  eV,  and  the  broad  dashes  represent  A Te  =  2  eV.  The  second  and 
third,  in  association  with  the  first  effect,  were  a  general  broadening  and  strengthen¬ 
ing  of  the  electric  field  and  a  variation  in  the  charge  distribution  around  the  shock 
front.  Figure  23  illustrates  the  variation  in  the  electric  field;  the  solid  line  represents 
the  constant  Te  =  2eV  case  and  the  dashed  line  represents  the  A Te  =  2eV  case. 
Figure  24  illustrates  the  net  charge  density,  which  was  approximated  via  Equation  2. 
In  this  case,  charge  neutrality  is  maintained  to  within  98.6  percent.  Again,  as  in  the 
previous  chapter,  the  magnitude  of  the  error  is  due  in  part  to  the  coarseness  of  the 
spatial  grid  used  in  the  numerical  solution.  For  all  values  of  A Te  and  6  investigated 
there  were  no  significant  effects  on  neutral  flow  variables,  shock  structure,  or  shock 
propagation,  which  is  expected  from  the  arguments  in  the  beginning  of  the  chapter. 
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Figure  23  Electric  field  variation  resulting  from  variable  electron  temperature. 
Constant  Te  =  2  eV  (solid  line)  and  A Te  =  2 eV  (dashes). 
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Figure  24  Approximate  net  charge  density  resulting  from  variable  electron  tem¬ 
perature.  A Te  =  2  eV  and  6  =  200. 
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4.1.3  Analysis.  The  increase  in  precursor  width  is  simply  the  result  of 
capturing  the  local  value  of  Te.  As  an  illustration,  recall  Avramenko’s  definition  of 
the  steady-state  precursor  width,  Equation  45: 

_  Vj  _  k_BTe 

CVin  mCVin 


The  precursor  width  varies  directly  with  electron  temperature  and  the  precursor 
exists  in  this  region  of  elevated  electron  temperature. 

The  profile  of  the  electric  field  is  affected  by  the  variable  electron  temperature. 
Recall  Equation  63: 


kB_ fTedrh  dTe. 
q  n*  d£  d£ 


When  Equation  64  is  subsituted  for  Te,  the  equation  becomes 

+  (67) 

q  rii  d£  oa  o  A 

The  term  on  the  left  is  generally  positive,  but  only  has  a  significant  value  near  £  -  0- 
The  term  on  the  right  is  positive  for  £  >  0,  and  negative  for  £  <  0,  which  creates 
the  negative  portion  of  the  electric  field.  For  f  <  0  and  £  0,  dpi/d£  and  the 

exponential  term  are  both  essentially  zero,  and  the  field  vanishes.  In  Figure  23  the 
broad  negative  region  behind  the  shock  draws  ions  away  from  the  shock  front,  further 
reducing  the  possibility  of  effective  ion-acoustic  wave  damping.  The  broad  positive 
field  region  imparts  energy  to  the  ions  and  extracts  energy  from  the  electrons,  but 
the  ion-neutral  coupling  at  this  fractional  ionization  is  not  strong  enough  to  affect 
the  neutral  population.  The  variation  in  ion  density  illustrated  in  Figure  22  shows 
the  influence  of  this  field,  particularly  in  the  partially  evacuated  region  just  behind 
the  shock.  The  total  potential  drop  over  the  region  is  about  -2  V ;  however,  since  net 
ion  production  is  still  restricted,  efficient  coupling  between  the  ions  and  neutrals  does 
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not  occur  and  the  neutral  flow  remains  unaffected.  Of  particular  note,  the  profile  of 
the  numerical  solution  to  the  electric  field  resembles  the  observed  profile  in  Sirghi’s 
research,  depicted  in  Figure  20,  including  a  depressed  region  of  the  field  for  £  <  0. 
This  correlation  between  observation  and  the  numerical  results  lends  credence  to  the 
incorporation  of  a  spatially-dependent  electron  temperature  in  Hilbun  s  code. 

4-2  Electron  Impact  Ionization 

4.2.1  Development.  In  a  glow  discharge,  the  plasma  is  generated  by  elec¬ 
tron  impact  ionization,  where  electrons  acquire  energy  from  the  applied  electric  field. 
Losses  generally  occur  through  diffusion  to  the  walls  of  the  tube.  As  stated  in  Chap¬ 
ter  I,  the  present  research  focuses  on  a  one-dimensional  approximation;  therefore, 
two-dimensional,  radial  loss  processes  cannot  be  considered.  The  production  and 
loss  of  plasma  are  assumed  to  balance  each  other.  This  section  investigates  the  ef¬ 
fects  of  allowing  additional  ionization  to  occur  at  the  shock  front.  Thus,  the  previous 
restriction  of  a  zero  net  production  mechanism  is  lifted.  The  growth  in  plasma  den¬ 
sity  is  only  restricted  by  the  energy  present  in  the  gas  and  by  the  assumed  discharge 
geometry.  Physical  mechanisms  to  mitigate  the  ionization  process  are  considered  in 
the  next  section. 

Even  though  Hilbun ’s  plasma  analysis  does  not  allow  for  additional  net  ioniza¬ 
tion,  he  concedes  that  there  is  some  experimental  evidence  that  additional  ionization 
is  present  at  the  shock  front  region:  “If  these  measurements  are  accurate,  then  it 
may  be  possible  for  the  fractional  ionization  in  the  shock  front  region  to  be  signif¬ 
icantly  greater  than  that  normally  present  in  the  quiescent  plasma  in  front  of  the 
shock.  In  such  a  case,  energy  transfer  from  the  ion-acoustic  wave  to  the  neutral 
shock  may  become  a  relevant  process,  resulting  in  a  perturbation  of  the  neutral  flow 
(22:98).”  The  research  to  which  he  refers  was  conducted  by  Gorshkov,  et  al.,  (19), 
and  Chutov  (13). 


49 


From  work  Chutov  had  completed  in  the  1970s,  he  discovered  that  “a  shock 
wave  in  a  gas-discharge  plasma  changes  the  degree  of  its  ionization  at  any  [shock] 
intensity.”  It  has  been  shown  that  the  presence  of  the  electronic  double  layer  due  to 
the  steep  gradients  at  the  shock  creates  a  potential  drop.  When  the  direction  of  the 
field  due  to  the  double  layer  coincides  with  that  of  the  electric  field  in  the  discharge, 
then  additional  energy  is  released  in  the  region  of  the  potential  drop,  which  leads 
to  additional  ionization  at  the  front.  Gorshkov  coined  this  effect  as  an  “electric 
detonation”  (13:506). 

Gorshkov  detected  additional  ionization  at  the  shock  front  in  decaying  glow 
discharges.  Shock  propagation  and  structure  were  measured  in  the  discharges  of 
a  variety  of  gases— Af,  N2,  and  CO2.  The  absence  of  an  ambient  electric  field 
during  passage  of  the  shock  is  similar  to  the  conditions  incorporated  into  the  model 
used  in  the  present  research.  Using  collimated  photomultipliers,  they  measured  the 
emissions  of  the  glow  discharge  plasma  and  observed  a  significant  drop  in  luminosity 
at  the  shock,  which  they  attributed  to  stimulated  ionization  of  excited  particles  at 
the  shock  front  (19:1141). 

The  present  research  will  assume  a  quasi-kinetic  approach  to  additional  ioniza¬ 
tion  at  the  shock  front.  The  number  of  ionization  events  per  unit  volume  per  unit 
timp  is  a  simple  result  of  kinetic  theory  assuming  a  Maxwellian  distribution  (e.g.  see 
(40:386)): 

dn  f°° 

zian  =  (-^)ion  =  nnne  J  ae(ve)vefe(ve)dve  =  nnneke  (68) 

where  ve  is  the  electron  thermal  velocity,  <7e(ve)  is  the  electron  impact  ionization 
cross  section,  and  /e(uc)  is  the  electron  velocity  probability  function.  The  integration 
extends  over  electron  thermal  speeds  whose  energies  exceed  the  ionization  potential; 
therefore  the  lower  limit  of  integration  is  Vk  =  \f^Ip/me.  Even  though  cre(ve) 
is  a  function  of  electron  energy,  for  the  purpose  of  approximation  it  is  assumed  to 
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be  roughly  constant.  For  all  cases,  <xe  was  given  an  order  of  magnitude  estimate 
of  10-21  m2.  The  Maxwellian  velocity  distribution,  fe(ve),  is  assumed  to  be  one¬ 
dimensional.  Therefore,  the  ionization  coefficient,  ke,  evaluates  to 
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The  ionization  coefficient  then  is  simply  a  function  of  electron  temperature.  In 
this  investigation,  the  electron  temperature  is  defined  by  Equation  64.  The  quantity 
IP/kB  in  the  Boltzmann  factor  can  be  considered  as  an  ionization  temperature,  TP. 
With  an  ionization  potential  of  15.8  eV  for  argon,  TP  evaluates  to  1.83  x  105  K-  This 
eliminates  consideration  of  ionization  by  ion-neutral  and  neutral-neutral  collisions. 
Since  these  particles  have  temperatures  on  the  order  of  300  K,  the  heavy  particle 
Boltzmann  factor  is  essentially  zero.  Free  electrons,  however,  yield  more  substantial 
ionization  coefficients.  For  A.Te/Teo  as  low  as  0.02,  the  ionization  coefficient  due  to 
electron  impact,  ke,  is  about  1.1  x  10-13 cm3/ s.  For  ambient  initial  conditions  of 
nn o  =  5.5  x  1018  cm-3  and  a  =  10-6,  the  ionization  rate  is  about  3.2  x  1012  cm-3  /is-1. 
The  average  computational  time  step  in  the  plasma  code  is  within  an  order  of  magni¬ 
tude  of  1  /xs;  therefore,  within  one  time  step,  a  60  percent  increase  in  the  ion  density 
at  the  shock  would  be  expected. 

A  similar  ionization  rate  as  that  above  could  be  derived  for  the  ambient  electron 
temperature  (2eV  or  23200  K);  therefore,  the  ionization  terms  must  be  balanced  or 
nullified  in  the  regions  away  from  the  shock.  This  balance  was  struck  by  using  the 
following  relationship: 


(  j_i  )net  —  Z-riet  —  U.nTle(A:e  ke o) 
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Therefore,  the  net  ionization  rate  is  simply  a  function  of  the  difference  in  the  expo¬ 
nential  terms.  Away  from  the  shock,  where  Te  =  Te o,  then  Z„et  =  0. 

Due  to  the  dependence  on  Equation  64,  the  electric  field  should  retain  the  basic 
characteristics  of  the  field  generated  in  the  previous  section  (Figure  23).  Since  the 
region  of  ionization  is  assumed  to  be  centered  on  the  shock  due  to  its  dependence 
on  Equation  64,  a  large  number  of  the  freed  electrons  would  tend  to  accumulate 
at  the  shock.  An  argument  was  presented  at  the  beginning  of  the  chapter  that 
promulgated  the  idea  that  the  newly  created  ions  would  generally  remain  at  the 
shock  front.  As  the  ion  density  builds  at  the  shock,  a  large  number  of  electrons 
should  be  expected  to  remain  in  the  vicinity  of  the  shock  as  well  due  to  the  mutual 
attraction  of  the  ions  and  electrons.  As  a  result  of  this  accumulation  of  charged 
particles,  an  exponential  increase  in  ion  density  in  time  is  expected  at  the  shock 
front,  as  determined  by  the  differential  equation,  Equation  70.  For  the  maximum 
electron  temperature  experienced,  Te  +  A Te,  Equation  70  can  be  integrated  to  find 
m(t )  at  the  shock  front,  assuming  quasi-neutrality,  n*  ~  ne,  applies: 
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Since  nn  is  several  orders  of  magnitude  larger  than  n,,  nn  can  be  considered  as  a 
constant  in  the  integration,  as  long  as  n{  only  comes  within  a  few  orders  of  magnitude 
of  nn: 
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rii(t)  =  niQ  exp(nn(ke  -  ke0)t) 


(73) 
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An  estimate  of  the  threshold  value  of  A Te  required  to  produce  a  neutral  shock 
acceleration  can  also  be  obtained.  In  the  unionized  Riemann  problem,  it  takes 
the  shock  approximately  tRiem(mn  =  0.97ms  to  propagate  from  x  —  0.4m  to  x  = 
0.9  m.  In  order  for  the  fractional  ionization  at  the  shock,  05 ,  to  reach  the  fractional 
ionization  required  for  effective  ion-neutral  momentum  coupling,  a{a,  within  tRiernann 
the  following  condition  must  be  met: 

'  Ip  s  <  Ip  Ip  \\  2/grjn (aja)  _ 

exp(  kBTe  V  rne  ^  BX^kBTe  kBT, *  ae  nntRiemann 

where  the  ionization  coefficients  have  been  expanded.  The  condition  yields  AT^/T^o 
0.01,  or  A Te  =  232  K,  as  a  lower  limit  to  observe  a  neutral  shock  acceleration  within 
the  time  hunt,  t Riemann  ~~  0.97  ms. 

Prom  Equation  73,  it  is  possible  to  estimate  the  time  for  as  to  reach  aia. 
Dividing  both  sides  by  nn 0  and  solving  for  t  yields: 

ln(a,./a)  (74) 

nn{ke  -  ke 0) 

For  ATe  =  300  K  and  a  —  10-6,  the  coefficients  evaluate  to  ke  =  9.76  x  10_20m3  /s 
and  fceo  =  8.82  x  10-20  m3/s,  which  gives  tia  ~  0.76  ms. 

Electron  impact  ionization  was  incorporated  into  Hilbun’s  one-dimensional 
plasma  code  in  much  the  way  that  the  variable  electron  temperature  was.  The  ion¬ 
ization  process  is  dependent  upon  the  electron  thermal  layer  at  the  shock;  therefore, 
the  variable  electron  temperature,  Equation  64,  was  used.  Incorporating  additional 
ionization  in  Hilbun’s  code  requires  the  modified  electric  field  source  term,  Equation 
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66,  and  the  introduction  of  a  new  ionization  source  vector: 

ZnetTTln 
Znet^vMn 

,  -Z^,(\v*  + 

04  = 
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Znet'W'iVi 

z^mQV?  + 

4.2.2  Results  and  Analysis.  For  a  range  of  the  parameters  ATe/Te0  and  6, 
a  number  of  test  cases  were  completed  for  the  initial  conditions  defined  previously. 
Three  trends  were  observed  in  neutral  shock  parameters,  as  exhibited  in  Figure  25. 
In  the  region  indicated  by  the  diamonds,  there  were  negligible  effects  on  neutral 
shock  parameters;  however,  the  electric  field  was  strengthened  and  broadened.  To¬ 
tal  electric  potentials  grew  to  —1.7  V,  indicative  of  the  influence  of  the  variation 
in  electron  temperature  at  the  shock.  In  the  parameter  zone  demarcated  by  the 
squares,  steady  accelerations  of  the  neutral  shocks  were  observed  with  an  increase 
in  electric  field  strength  and  width.  Shock  accelerations  were  observed  for  values 
of  ATe/Teo  just  greater  than  the  threshold  determined  above.  Incidentally,  shock 
acceleration  generally  occurred  when  the  potential  reached  — 6  V.  There  were  also 
modifications  of  the  neutral  shock  profile  that  could  be  interpreted  as  a  weakening  of 
shock  strength.  Perhaps  the  most  intriguing  result  was  the  appearance  of  a  signif¬ 
icant  neutral  precursor,  in  the  regime  marked  by  the  triangles.  For  convenience,  a 
significant  precursor  is  defined  here  as  a  10  percent  rise  in  upstream  pressure  in  front 
of  the  shock.  The  strengths  and  spatial  extents  of  the  electron  thermal  layer  em¬ 
ployed  in  these  test  runs  are  based  on  empirical  results  rather  than  a  solution  of  the 
electron  continuity,  momentum,  and  energy  equations.  Temperature  increases  and 
gradients  that  were  too  large  led  to  nonphysical  conditions  and  unstable  solutions. 
The  region  indicated  by  the  crosses  represents  these  cases. 
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Figure  25  Results  of  unrestricted  electron  impact  ionization  at  the  shock  front. 

The  parameter  space  is  represented  by  the  horizontal  axis  in  terms  of 
ATe/Teo,  and  the  vertical  axis  in  terms  of  8,  measured  in  neutral-neutral 
mean  free  paths.  Lower  parameter  values  yielded  strengthened  and 
broadened  electric  fields.  Higher  values  produced  shock  accelerations 
and  neutral  precursors. 
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Figure  26  Acceleration  of  shock  due  to  additional  ionization. 

Acceleration  of  the  neutral  shock  was  generally  observed  for  A Te  >  300  K  (for 
Te0  =  2  eV) ,  just  above  the  threshold  value  of  A Te  derived  above.  The  parameters  of 
the  case  discussed  below  are  A Te  =  300  K  and  <5  =  300.  The  acceleration  begins  at 
approximately  t  —  0.45  ms  after  rupture  of  the  diaphragm,  as  illustrated  in  Figure  26, 
which  is  about  half  the  value  of  tia  predicted  by  Equation  74  for  A Te  =  300  K-  That 
prediction,  however,  did  not  take  the  width  of  the  ionization  region,  6,  into  account. 
If  a  linear  velocity  relationship  is  assumed  over  the  time  period  of  acceleration,  the 
shock  is  seen  to  accelerate  at  approximately  1.4  x  105  m/s2. 

The  neutral  shock  pressure  and  density  profiles  shown  in  Figures  27  and  28, 
respectively,  do  not  match  the  profiles  found  in  literature  ((20),  (27),  (31),  and  (30)). 
Some  distance  behind  the  shock,  there  is  a  weakening  in  the  pressure.  Whereas  the 
pressure  ratio  of  the  control  case  is  2.94,  the  lowest  pressure  ratio  experienced  in  this 
region  is  P2/P1  =  2.22,  which  is  a  25  percent  reduction.  On  the  other  hand,  as  seen 
in  the  figures,  there  is  also  a  strengthening  of  the  shock.  For  the  case  illustrated, 
Pi/ Pi  grows  from  2.94  to  4.27.  From  Equation  16,  a  pressure  ratio  of  this  magnitude 
corresponds  to  a  jump  in  Mach  number  from  1.60  to  1.90,  or  a  growth  in  shock  speed 
from  516  m  /  s  to  613  m  /  s.  The  shock  speed,  however,  accelerates  to  only  590  m  /  s, 
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Figure  27  Modification  of  neutral  pressure  profile  for  an  accelerating  shock  (solid) 
due  to  additional  ionization,  6  =  300  and  A Te  =  300  K-  Control  case  is 
dashed.  Normalized  to  upstream  pressure.  Weakening  and  broadening 
of  the  shock  was  not  observed. 

as  seen  in  Figure  26.  It  falls  short  of  the  anticipated  613  m/s  possibly  because  the 
shock  no  longer  resembles  a  Riemann  shock,  meaning  that  it  no  longer  has  the  level 
pressure  reservoir.  This  modification  of  the  neutral  shock  structure  is  the  result 
of  the  transport  of  neutrals  from  Region  2  to  the  shock  front  due  to  ion-neutral 
coupling,  as  made  evident  by  the  shift  in  neutral  density  in  Figure  28. 

Equation  73  predicted  the  exponential  growth  of  ion  density  with  time.  When 
the  shock  had  reached  x  =  0.9  m,  the  ion  density  in  the  vicinity  of  the  shock  had 
increased  by  four  orders  of  magnitude,  as  illustrated  in  Figure  29.  From  Figure  26, 
the  shock  experienced  an  appreciable  acceleration  at  t  =  0.45  ms,  which  corresponds 
to  a  local  fractional  ionization  at  the  shock  front  of  approximately  1  x  10-3,  as  seen  in 
Figure  30.  This  is  in  very  good  agreement  with  the  threshold  fractional  ionization, 
cti a,  described  in  the  beginning  of  the  chapter. 

It  is  possible  to  approximate  the  shock  speed  as  a  function  of  time.  In  Figure  30 
as  is  shown  as  a  function  of  time.  In  Figure  17  the  shock  speed  is  shown  as  function 
of  as-  If  the  gas  was  fully  ionized,  the  fractional  ionization  would  be  unity  and 
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Figure  28  Modification  of  neutral  density  profile  for  an  accelerating  shock  (solid) 
due  to  additional  ionization,  6  =  300  and  A Te  =  300  K-  Control  case  is 
dashed.  Normalized  to  upstream  density.  Weakening  and  broadening 
of  the  shock  was  not  observed. 


Distance  (cm) 

Figure  29  Ion  density  due  to  additional  ionization  at  the  shock,  6  =  300  and 
A Te  —  300  K-  Density  is  normalized  to  upstream  ion  density.  In 
the  vicinity  of  the  shock,  the  local  value  of  a  grows  to  approximately 
2  x  10-2,  which  results  in  a  bidirectional  pressure  gradient  in  Region  2. 
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Figure  30  Growth  of  local  fractional  ionization  at  the  shock  front  due  to  additional 
ionization. 

the  shock  would  propagate  at  the  ambient  ion-acoustic  velocity,  which  is  2196m /s 
according  to  Equation  23  for  test  conditions.  Incorporating  this  data  point  into  the 
data  shown  in  Figure  17,  the  data  can  be  approximated  by  the  following  curve: 

c  =  co(l  +  3.28c*a676)  (76) 

where  Co  is  the  shock  velocity  for  the  unionized  Riemann  solution,  namely,  in  this 
case,  513  m  /  s.  The  data  in  Figure  30  can  be  approximated  by 

as  =  0.02718(— )4'323  (77) 

lms 

Combining  these  two  relationships  yields 

c  =  co(l  +  0.287(-^— )2,92)  (78) 

v  lms 

This  curve  is  plotted  along  with  the  shock  propagation  data  in  Figure  31.  Although 
the  curve  follows  the  general  trend  of  the  numerical  solution,  the  curve  predicts  a 
much  higher  velocity  than  the  results  of  the  numerical  solution.  The  disparity  could 
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Figure  31  Acceleration  of  shock  due  to  additional  ionization  compared  to  predic¬ 
tion. 

be  the  result  of  estimating  two  exponential  functions  from  a  limited  number  of  data 
points,  or  the  result  of  the  assumption  that  the  values  of  as  are  always  exactly  at 
the  neutral  shock  front,  which  is  not  always  the  case  as  seen  in  Figure  29. 

In  spite  of  the  exponential  increase  in  ion  density  evident  in  Figure  29,  the 
ion  density  profile  no  longer  resembles  that  of  a  Riemann  shock.  In  addition  to 
the  pressure  gradient  associated  with  the  shock,  there  is  also  a  negative  pressure 
gradient  that  causes  ions  to  move  against  the  direction  of  shock  propagation  and 
into  Regions  2  and  3.  This  also  results  in  an  extended  negative  electric  field  region 
behind  the  shock,  which  extends  to  the  contact  discontinuity,  as  seen  in  Figure  32. 
The  sharp  negative  spike  in  the  field  may  be  an  artifact  of  the  numerical  solution, 
since  it  was  necessary  to  inhibited  electric  field  effects  in  the  code  behind  the  contact 
discontinuity.  Even  though  the  concern  of  the  present  research  is  with  the  field  in 
the  vicinity  of  the  shock,  the  total  potential  is  obtained  by  integrating  over  the 
entire  negative  and  positive  field  regions,  just  as  in  Section  4.1.3.  The  potential 
drop  delivered  by  the  positive  field  is  approximately  —19  V  and  the  potential  drop 
of  the  negative  portion  is  +11 V;  therefore,  giving  a  total  potential  drop  of  — 8  V. 
The  evolution  of  the  potential  as  the  shock  propagates  and  the  ion  density  builds 
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Figure  32  Growth  and  broadening  of  electric  field  due  to  additional  ionization  at 
the  shock,  8  =  300  and  A Te  =  300  K.  Total  potential  due  to  this  field 
is  approximately  -8  V. 

is  illustrated  in  Figure  33.  Incidentally,  shocks  generally  began  to  accelerate  when 
they  experienced  a  potential  drop  of  —6  V*  This  high  potential  cannot  be  the  cause 
of  the  acceleration  as  shown  above;  however,  the  potential  enhances  the  ion-acoustic 
waves  by  way  of  accelerating  the  charged  particles.  The  electric  field  strength  was 
observed  to  increase  with  additional  ionization  over  the  control  case  by  a  factor  of 
2.5,  and  the  width  increased  by  a  factor  of  5  (compare  to  Figure  15). 

Recall  from  Section  1.3  that  Saeks  and  Kunhardt  (35)  proposed  that  the  equa¬ 
tion  of  state  of  a  weakly  ionized  gas  should  include  electrostatic  pressure  terms.  The 
electrostatic  pressure  is  given  by 

Pe  =  \e0  E2  (79) 

where  the  permittivity  of  free  space,  eo>  equals  8.85  x  10~12  C2  /  N  m2-  Given  the 
maximum  field  strength  in  Figure  32  of  approximately  5  V  /  cm,  the  maximum  elec¬ 
trostatic  pressure  is  about  1.5  x  10~4  torr.  Compared  to  the  upstream  ambient 
neutral  thermal  pressure  of  30  torr,  the  electrostatic  pressure  is  negligible.  Now 
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Figure  33  Growth  of  total  potential  as  a  function  of  time  due  to  additional  ioniza¬ 
tion  at  the  shock,  6  =  300  and  A Te  =  300  K- 

compare  the  the  electrostatic  pressure  gradient  force  on  the  charged  component  to 
the  pressure  gradient  force  of  the  shock.  The  electrostatic  pressure  gradient  is  given 

by 

(IPe  /'on'\ 

T—** I  m 

which  yields  a  pressure  gradient  force  per  unit  mass  of 


€qE  dE 


(81) 


In  the  vicinity  of  the  shock,  the  average  field  strength  is  3  V/ cm,  the  average  field 
gradient  of  1 V  /  cm2 ,  and  the  ion  mass  density  is  on  the  order  of  10~8  kg  /  cm3 .  This 
yields  an  average  force  per  mass  due  to  the  electric  field  on  the  ions  is  on  the  order 
of  10-2  N  /  kg.  However,  the  pressure  gradient  force  per  mass  due  to  the  shock  on 
the  ions  is  on  the  order  of  105  N  /  kg.  The  pressure  gradient  force  due  to  the  electric 
field  is  insignificant  when  compared  to  that  of  the  thermal  and  dynamic  pressure 
gradients  of  the  gas  as  a  whole  and  should  not  be  given  any  consideration  in  the 
equation  of  state  of  weakly  ionized  gases. 
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Figure  34  Net  charge  density.  Growth  and  broadening  of  electronic  double  layer 
due  to  additional  ionization  at  the  shock,  6  =  300  and  A Te  =  300  K- 

The  net  charge  density  in  the  vicinity  of  the  shock  is  displayed  in  Figure  34. 
Charge  neutrality  has  been  maintained  to  99.2  percent.  The  long  negative  tail  of 
the  electric  field  serves  a  conduit  to  draw  ions  away  from  the  highly  concentrated 
region  at  the  shock  front  and  into  the  remainder  of  Region  2,  thus  expanding  the 
EDL  (compare  to  Figures  12  and  16).  This  is  visible  in  Figure  35,  in  which  the  ion 
velocity  downstream  from  the  shock  is  reduced  to  75  percent  of  that  of  the  control 
case. 

Significant  neutral  precursors  were  observed  for  cases  in  which  ATe  >  500  K. 
For  the  particular  case  discussed  below,  A Te  =  600  K  and  6  =  100.  The  precur¬ 
sor  is  clearly  visible  in  Figures  36  and  37,  with  a  width  of  about  a  centimeter  and 
a  maximum  pressure  of  approximately  1.1PX.  Increased  density  in  the  precursor 
indicates  mass  transport  from  Region  2  to  Region  1,  which  must  be  the  result  of 
ion-neutral  momentum  coupling.  These  precursors  proved  to  be  transient  phenom¬ 
ena.  The  precursor  in  Figure  36  occurred  at  t  =  0.43  ms,  where  it  had  reached  its 
maximum  intensity,  well  before  the  shock  had  reached  x  =  0.9  m.  By  the  time  the 
shock  had  reached  x  =  0.9  m,  the  precursor  had  diminished  below  the  criteria  estab¬ 
lished  at  the  beginning  of  this  section.  In  Hilbun’s  plasma  code,  neutral  pressure, 
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Figure  35  Ion  velocity  due  to  additional  ionization  at  the  shock,  6  =  300  and 

ATe  =  300 K-  Control  case  shown  for  comparison  (dashed).  Velocities 
are  normalized  to  upstream  ion-neutral  thermal  velocity,  v  =  250 m/s. 
Ions  are  slowed  by  the  large  negative  pressure  gradient  in  the  region  of 
ionization. 

Pn,  is  normalized  by  the  ambient  upstream  pressure,  Pn o,  and  all  calculations  use 
pn  =  Pn/Pn0.  Shock  detection  is  facilitated  by  starting  from  x  =  1  m  and  searching 
down  the  shock  tube  until  the  first  pressure  rise,  P„  >  Peru  ,  is  detected,  where  Per  it 
is  arbitrary  value  just  above  unity.  Precursors  were  detected  by  setting  Parit  =  1-1- 

In  conclusion,  unrestricted  electron  impact  ionization  sustained  by  a  localized 
increase  in  electron  temperature  at  the  shock  front  can  modify  neutral  shock  struc¬ 
ture  arid  propagation.  Additional  ionization  is  effective  in  enhancing  ion-acoustic 
wave  damping  by  raising  the  local  fractional  ionization  at  the  shock  front  to  a  level 
conducive  to  ion-neutral  momentum  coupling,  such  that,  at  the  shock,  the  medium 
becomes  a  partially  ionized  gas,  where  a  >  10"3.  However,  research  indicates 
that  the  shock  acceleration  in  weakly  ionized  gases  does  not  continue  unhampered; 
therefore,  there  must  be  some  abating  mechanism  present. 
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Figure  36  Neutral  Precursor.  Neutral  pressure  profile  at  the  shock  front  for  the 
parameters  A Te  =  600  K  and  6  =  100.  Precursor  is  clearly  visible 
against  the  control  case  (dashed). 


Distance  (cm) 

Figure  37  Neutral  Precursor.  Neutral  density  profile  at  the  shock  front  for  the 
parameters  A Te  =  600  K  and  S  =  100.  Precursor  is  clearly  visible 
against  the  control  case  (dashed).  Increased  density  in  the  precursor 
indicates  mass  transport  from  Region  2  to  Region  1. 
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Figure  38  Shock  speed  as  a  function  of  time  for  a  shock  in  N2  (dotted)  in  Ar 
(solid).  Shock  enters  plasma  region  at  t  =  0  (Ref:  (12)). 

4.3  Volumetric  Loss  of  Plasma  at  the  Shock 

4.3.1  Development.  Revisiting  the  shock  velocity  depicted  in  Figure  26, 
all  accelerating  shocks  observed  in  the  current  research  continued  to  accelerate  until 
they  reached  the  end  of  the  shock  tube.  Research,  however,  indicates  that  shocks 
level  off  in  velocity  after  their  initial  acceleration  ((12),  (13),  (27)).  Figure  38  displays 
the  experimental  results  obtained  by  Chutov,  et  al.,(12)  in  which  the  acceleration 
of  the  shock  leveled  off  when  it  entered  the  a  region  of  weakly  ionized  argon.  For 
comparison,  their  results  for  weakly  ionized  N2  are  included.  Chutov  attributes  this 
non-monotonic  acceleration  in  N2  to  vibrational-translational  interaction,  which  is 
not  present  in  argon.  It  is  apparent  then  that  there  is  some  mechanism  at  work  to 
mitigate  the  unabated  acceleration  of  the  shock.  This  section  considers  volumetric 
loss  of  plasma  due  to  ion-electron  recombination  as  a  candidate  mechanism. 

In  the  absence  of  applied  electric  fields  and  negative  ions,  the  ion-electron 
recombination  equation  is  given  by  Raizer  (34:60)  as 

Zrecamb  =  )recawb  =  -rtiUePe  (82) 


66 


where  fie  is  the  dissociative  recombination  coefficient,  which  is  generally  on  the  order 
of  10— 7  cm3  /  s  for  diatomic  gaseous  media.  The  net  rate  of  ionization  events  per 
unit  volume  per  unit  time,  Z^t ,  is  the  sum  of  Equations  68  and  82: 

Znet  =  Zion  ~  Z recomb  =  ne{keTln  —  PeTli)  (83) 


Note  that  the  ion  loss  term  now  has  an  nj  dependence,  rather  than  an  n*  depen¬ 
dence  as  in  the  unrestricted  ionization  cases  above.  In  order  to  restrict  net  impact 
ionization  to  vicinity  of  the  neutral  shock  front,  Znet  must  be  positive  in  the  vicinity 
of  the  neutral  shock  and  effectively  zero  for  £  0  and  £  0,  where  tiiq  otnn o- 

For  this  condition  to  be  met 
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Raizer  states  that  for  the  dissociative  recombination  process  of  argon  ions, 


Ar\  +  e  — »  Ar  +  Ar* 


Pe  varies  as  Te_1/2  for  temperatures  from  room  temperature  to  several  thousand 
degrees  Kelvin,  and  varies  as  Te-3/2  for  even  higher  temperatures.  As  an  approxi¬ 
mation,  the  present  research  will  assume  that  (3e  varies  as  Te  1,  so  that 
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(Te)  =  Pe 0 


Te 0 
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(85) 


Therefore,  where  Te  —  Te o,  then  (3e  —  Pe o  —  keo/a,  and  the  net  ionization  is  zero. 
Electron  impact  ionization  mitigated  by  ion-electron  recombination  was  also  imple¬ 
mented  in  Hilbun’s  plasma  code  by  incorporating  Equations  64,  66,  75,  and  83,  where 
quasi-neutrality,  n;  ^  ne,  is  assumed. 
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Figure  39  Ion  density  profile  due  to  additional  ionization  mitigated  by  and  re¬ 
combination  at  the  shock,  8  =  300  and  A Te  =  8000  K-  Density  is 
normalized  to  upstream  ion  density.  In  the  vicinity  of  the  shock,  the 
local  value  of  a  grew  to  approximately  2.7  x  10-5  and  remained  there. 

4.3.2  Results  and  Analysis.  Calculations  with  input  parameters  in  the 
range  of  A Te  <  11, 000  K  and  8  <  300  were  in  general  numerically  stable.  Dramatic 
effects  on  ion  and  neutral  flow  variables  were  not  observed,  neither  were  significant 
accelerations  of  the  shock  front.  The  parameters  of  the  case  discussed  below  are 
ATe  =  8000  K  and  8  =  300.  In  Figure  39,  the  ion  density  at  the  shock  front  increased 
such  that  as  =  2.5  x  10-5,  where  it  remained  for  the  duration  of  its  propagation. 
This  value  is  greatly  diminished  from  the  unrestricted  ionization  values  due  to  the 
n2  recombination  factor.  This  value  of  as  does  not  approach  the  1  x  10-3  required 
for  effective  ion-acoustic  wave  damping.  Consequently,  there  were  no  variations  in 
the  neutral  shock  pressure  profile.  The  neutral  shock  velocity  was  seen  to  increase 
only  very  slightly,  if  at  all.  The  general  trend  in  the  shock  velocity  was  a  change  of 
2  m  /  s  over  0.5  ms;  however,  the  uncertainty  in  the  shock  velocity  calculations  was 
also  ±2  m  /  s.  Acceleration  of  the  shock  in  this  case  is  slight  at  best. 
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Figure  40  Critical  fractional  ionization  at  the  shock  front  as  a  function  of  the 
increase  in  electron  temperature. 

The  growth  of  ion  density  at  the  shock  is  limited  by  the  ion-electron  recombi¬ 
nation.  Recall  Equation  83: 


Znet  —  Z ^0^  Zrecomt)  —  Tle{keTln  PeP'i ) 


Ion  density  grows  until  Znet  =  0,  at  which  point  rii  —  (ke//3e)nn,  if  it  is  assumed  that 
nn  is  roughly  constant  during  this  process.  A  cutoff  as  can  be  determined  from 
Equations  69,  71,  and  85: 


aSc  =  -T- 


ke 

We 


=  a(l  +  ^)exp(  Ip 

e0 


(86) 


This  function  is  plotted  in  Figure  40.  Clearly,  a  large  increase  in  electron  tem¬ 
perature  is  required  to  reach  asc  —  aia ■  These  increases  in  electron  temperature, 
however,  caused  numerical  instabilities. 

The  electric  field,  depicted  in  Figure  41,  exhibited  the  general  profile  of  the 
variable  electron  temperature  case  shown  in  Figure  23.  This  is  to  be  expected  since 
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Figure  41  Electric  field  in  vicinity  of  shock  due  to  additional  ionization  mitigated 
by  recombination  at  the  shock,  6  =  300  and  A Te  =  8000  K- 


the  growth  of  ion  density  at  the  shock  is  limited,  yet  the  field  is  still  affected  by  the 
variation  in  electron  temperature.  The  field  reaches  a  maximum  of  approximately 
4  V  /  cm  at  the  shock  boundary.  It  also  has  a  broad  positive  precursor  and  negative 
tail,  both  features  associated  with  the  electron  temperature  dependence  of  Equation 
67.  The  broad  positive  field  region  energizes  the  charged  components,  but  there 
is  insufficient  ion  density  to  affect  the  neutral  population.  The  broad  negative 
region  behind  the  shock  serves  as  a  conduit  to  draw  ions  away  from  the  shock  front, 
further  reducing  the  possibility  of  effective  ion-neutral  momentum  coupling.  Here, 
the  influence  is  even  more  profound  as  seen  in  Figure  41.  The  potential  drop  of  the 
positive  field  is  -8  V,  and  the  potential  in  Region  2  is  +6.5  V.  In  spite  of  the  broad 
extent  of  the  field,  the  total  potential  drop  delivered  by  the  field  is  only  -1.5  V. 

The  net  charge  density  also  exhibited  some  interesting  characteristics,  as  seen 
in  Figure  42.  In  the  absence  of  the  spike  in  the  electric  field,  the  lobes  would  tend  to 
trap  electrons  at  the  shock  front.  With  the  spike  present,  electrons  are  trapped  in  a 
small  region  just  in  front  of  the  shock  as  seen  in  Figure  42,  in  contrast  to  the  broad 
negative  charge  regions  of  the  previous  cases.  The  field  structure  also  establishes 
a  double  EDL,  in  which  there  are  two  regions  of  net  positive  charge  near  the  shock 
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Figure  42  Net  charge  density  in  vicinity  of  shock  due  to  additional  ionization  mit¬ 
igated  by  and  recombination  at  the  shock,  8  =  300  and  A Te  —  8000  K- 

front.  In  spite  of  the  general  negative  appearance  of  the  net  charge  profile,  charge 
neutrality  was  maintained  to  within  97.7  percent. 

Due  to  the  lack  of  effects  on  the  neutral  shock  structure  and  propagation,  it 
is  likely  that  either  ion-electron  recombination  occurs  on  a  smaller  scale  than  that 
which  occurred  in  this  investigation,  or  that  there  are  other  mechanisms  at  work  in 
the  plasma.  Some  of  these  mechanisms  may  be  hidden  in  the  limitations  imposed 
upon  the  model,  which  were  discussed  in  Section  1.3.  These  may  include  energy 
losses  to  electronic  excitation,  diffusion  losses,  and  modifications  to  electron  energy 
distribution  functions. 
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V.  Conclusions  and  Recommendations 

5.1  Conclusions 

In  Chapter  II,  several  phenomena  associated  with  shock  propagation  in  weakly 
ionized  gases  were  described.  Four  of  those  were  investigated  in  this  research: 

•  Shock  waves  have  an  anomalously  high  propagation  velocity. 

•  The  shock  front  is  significantly  broadened. 

•  A  precursor  exists  ahead  of  the  shock  wave. 

•  Shock  strength  is  reduced. 

The  introduction  of  a  spatially-dependent  electron  temperature  profile  greatly 
affected  the  charged  component  for  the  values  of  A Te  and  6  investigated.  Three 
effects  on  the  charged  component  flow  parameters  and  electric  field  structure  resulted 
from  incorporating  this  temperature  profile.  The  first  was  an  increase  in  the  charged 
component  precursor  width,  £0.  This  was  shown  to  be  the  result  of  capturing 
the  local  value  of  Te,  since  £0  varies  directly  with  Te.  The  second  was  a  general 
broadening  and  strengthening  of  the  electric  field.  This,  however,  included  a  negative 
field  region  that  prevented  the  total  potential  drop  from  growing  very  large,  resulting 
in  a  drop  of  only  -2  V,  whereas  the  potential  drop  of  the  unionized  Riemann  value 
was  -1.5  V.  The  third  was  a  variation  in  the  net  charge  distribution  around  the 
shock  front  as  a  result  of  the  electric  field.  All  of  these  effects  can  be  explained 
by  the  electric  field  approximation’s  dependence  on  the  spatially-dependent  electron 
temperature  (Equation  67).  There  were  no  discernible  effects  on  the  neutral  gas, 
which  was  expected  since  the  ion  density  never  approaches  the  level  required  for 
efficient  ion-neutral  momentum  coupling. 

Although  the  spatially-dependent  electron  temperature  alone  bore  little  influ¬ 
ence  on  the  neutral  flow,  when  coupled  with  additional  electron  impact  ionization 
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at  the  shock  front,  it  yielded  striking  results.  For  the  range  of  the  parameters  A Te 
and  6  investigated,  three  trends  were  observed  in  neutral  shock  parameters.  Steady 
accelerations  of  the  neutral  shocks  were  observed.  There  were  also  modifications  of 
the  neutral  shock  profile.  Perhaps  the  most  intriguing  result  was  the  appearance 
of  a  significant  neutral  precursor.  These  results  were  triggered  by  an  increase  in 
ion  density  at  the  shock  to  the  level  required  for  sufficient  ion-neutral  momentum 
coupling.  It  was  determined  that  in  order  for  ion-acoustic  wave  damping  to  affect 
the  neutral  flow,  the  local  value  of  a  at  the  shock  front  must  increase  to  1  x  10-3. 
Given  this  condition,  a  time,  Ua,  and  a  minimum  required  A Te  were  derived  for 
expected  accelerations  of  the  shock.  Accelerations  were  observed  for  values  of  A Te 
just  greater  than  the  threshold  and  occurred  at  times  that  were  half  of  tia,  as  seen 
in  Figure  26. 

In  Region  2  of  the  shock  tube,  the  neutral  pressure  profile  deviated  significantly 
from  the  Riemann  solution.  A  general  weakening  of  the  cross-shock  pressure  ratio, 
P2/P1,  was  observed.  Although,  locally  at  the  shock  front,  there  was  also  an  increase 
in  P2/Pi.  This  modification  of  the  neutral  shock  structure  is  the  result  of  the 
transport  of  neutrals  from  Region  2  to  the  shock  front  due  to  ion-neutral  coupling, 
as  made  evident  by  the  shift  in  neutral  density  in  Figure  28. 

Significant  neutral  precursors  were  observed  for  higher  values  of  ATe.  The 
precursor  is  clearly  visible  in  Figure  36,  with  a  width  of  about  a  centimeter  and  a 
maximum  pressure  of  approximately  1.1  Pi.  These  precursors  proved  to  be  transient 
phenomena,  in  which  most  precursors  remained  above  Pa-u  for  only  0.03  ms. 

The  striking  modifications  to  the  neutral  shock  structure  and  propagation  due 
to  unrestricted  ionization  were  mitigated  by  ion-electron  recombination.  Due  to  the 
nf  recombination  factor,  the  local  value  of  a  did  not  approach  the  1  x  10-3  required 
for  effective  ion-acoustic  wave  damping.  The  electric  field  exhibited  the  structure 
expected  from  its  dependence  on  the  variable  electron  temperature.  The  positive 
portion  of  the  field  exhibited  two  maxima.  The  smaller  maximum  is  balanced  by  the 
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negative  portion  of  the  field  behind  the  shock,  which  serves  as  a  conduit  to  draw  ions 
away  from  the  shock  front,  further  reducing  the  possibility  of  effective  ion-acoustic 
wave  damping. 

In  conclusion,  it  has  been  shown  that  unrestricted  electron  impact  ionization 
sustained  by  a  localized  increase  in  electron  temperature  at  the  shock  front  can 
modify  neutral  shock  structure  and  propagation.  Additional  ionization  is  effective 
in  enhancing  ion-acoustic  wave  damping  by  raising  the  local  fractional  ionization 
at  the  shock  front  to  a  level  conducive  to  ion-neutral  momentum  coupling,  such 
that,  at  the  shock,  the  medium  becomes  a  partially  ionized  gas,  where  a  >  10  3. 
However,  research  indicates  that  the  shock  acceleration  in  weakly  ionized  gases  does 
not  continue  unhampered;  therefore,  there  must  be  some  abating  mechanism  present 
in  weakly  ionized  gases. 

5.2  Recommendations  for  future  study 

The  present  research  lifted  two  of  the  restrictions  imposed  upon  the  previous 
analytical  and  numerical  treatments.  There  are  many  remaining  to  be  conquered. 
Here  they  are  discussed  in  order  of  feasibility.  Many  of  the  following  physical 
considerations  could  feasibly  be  incorporated  into  Hilbun’s  codes.  The  numerical 
solutions  of  the  present  research  essentially  rode  on  the  shoulders  of  Hilbun’s  work; 
possible  improvements  to  that  work  are  also  outlined  below. 

Shocks  are  multi-dimensional  phenomena.  The  problems  studied  in  the  present 
research  are  limited  to  one  spatial  dimension.  In  his  research  on  the  effects  of 
thermal  inhomogeneities  on  shock  propagation,  Hilbun  developed  a  two-dimensional 
fluid  code  (22:28).  Plasma  effects  could  be  expanded  into  two  dimensions  and 
incorporated  into  this  program. 

Although  this  research  effort  focused  on  shock  tube  problems,  much  of  current 
research  is  conducted  in  glow  discharge  tubes.  In  a  glow  discharge  tube,  an  electric 
field  is  applied  to  the  gaseous  medium,  which  causes  electronic  excitation  and  ioniza- 
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tion.  The  ionization  is  balanced  by  diffusion  losses  to  the  tube  walls.  The  constant 
baseline  electric  field  of  the  positive  column  was  not  a  consideration  in  the  present 
research,  nor  were  radial  losses.  After  modifying  Hilbun’s  one-dimensional  plasma 
code  for  two  dimensions,  glow  discharge  tube  problems  could  be  more  accurately 
modelled. 

Any  practical  investigation  of  plasma  aerodynamics  must  include  diatomic 
gases.  With  these  species  come  a  host  of  considerations:  the  energies  associated 
with  rotation,  vibration,  and  dissociation.  In  his  study  of  post-shock  energy  addition 
due  to  vibrational  relaxation,  Hilbun  developed  a  one-dimensional,  time-dependent, 
single  fluid  code  (22:103)  that  takes  vibrational  energy  into  account.  With  great 
meticulousness,  his  plasma  and  vibrational  codes  could  be  melded  into  one.  Plasma 
effects  could  be  introduced  into  this  code.  Along  these  lines,  electronic  excitation 
should  be  considered  as  well  since  much  energy  can  be  absorbed  into  excited  states, 
thus  reducing  the  possibility  of  additional  ionization  occurring  at  the  shock  front. 

Several  assumptions  were  made  in  modelling  the  electron  impact  ionization 
process.  One  of  these  assumptions  was  a  steady  electron  energy  distribution  func¬ 
tion.  Another  assumption  was  a  constant  ionization  cross  section.  As  free  electrons 
lose  energy  to  the  ionization  process,  a  shift  in  the  electron  energy  distribution  func¬ 
tion  occurs.  Since  additional  ionization  was  shown  to  affect  neutral  shock  structure 
and  propagation  characteristics,  a  purely  kinetic  approach  to  the  ionization  process 
is  warranted. 

Charge  density  gradients,  according  to  the  steady-state  and  time-dependent, 
two-fluid  approximations,  vary  rapidly  just  behind  the  shock.  A  fixed,  laboratory 
spatial  scale  may  not  accurately  capture  these  gradients.  A  variable  spatial  scale 
could  be  implemented  in  the  vicinity  of  the  shock.  For  this,  the  fluid  equations 
must  be  transformed  from  the  laboratory-fixed  frame  to  the  shock-fixed  frame.  Ac¬ 
cording  to  Anderson  (3:102),  the  differencing  of  the  FCT  algorithm  can  be  reversed 
to  transform  to  a  shock-fixed  frame.  The  real  challenge  is  that  shock  propagation 
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in  the  Riemann  problem  is  nonlinear  at  the  onset  of  propagation.  An  algorithm  to 
predict  the  location  of  £  =  0  would  need  to  be  incorporated  as  well. 

Finally,  the  two-fluid  approximation  still  stands.  This  restriction  prevents  ions 
and  electrons  from  moving  independently  of  each  other ,  which  is  necessary  in  order 
to  fully  characterize  the  double  layer  influence  on  shock  structure  and  propagation 
in  weakly  ionized  gases.  The  electric  field  is  an  approximation  based  upon  the 
assumption  of  steady-state  electron  momentum.  In  order  to  fully  characterize  the 
electric  field  and  associated  charge  separation,  Poisson’s  equation  must  be  solved. 
Appendix  B  outlines  the  development  of  the  fluid  equations  to  this  end;  however,  a 
numerical  scheme  that  can  mitigate  the  restrictive  time  step  is  still  required. 
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Appendix  A.  List  of  Symbols 

The  subscripts  j  and  k  refer  to  the  neutral,  ion,  and  electron  species  as  appropriate. 


a 

speed  of  sound 

c 

shock  speed 

cp 

specific  heat  (constant  pressure) 

cv 

specific  heat  (constant  volume) 

e 

elementary  charge 

E 

electric  field 

fe 

electron  velocity  distribution 

F 

flux  vector 

Ip 

ionization  potential 

k 

ion-acoustic  wave  number 

kg 

Boltzmann’s  constant 

ke 

electron  impact  ionization  coefficient 

m 

ion-neutral  mass 

me 

electron  mass 

M 

Mach  number 

Tlj 

species  number  density 

Ncfl 

Courant-Priedrichs-Levy  number 

Pi 

species  pressure 

Pi 

species  normalized  pressure 

Per  it 

normalized  shock  detection  pressure 

Pjk 

inter-species  momentum  transfer 

Q 

charge 

Q 

heat  addition 

Qjk 

inter-species  energy  transfer 

s 

entropy 
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si 

Euler  term  source  vector 

S~2 

species  coupling  source  vector 

S~3 

electric  field  source  vector 

si 

ionization  source  vector 

t 

time  coordinate 

tia 

time  to  reach  ion-acoustic  conditions 

Tj 

temperature 

TP 

ionization  temperature 

u 

fluid  velocity 

U 

vector  of  conserved  variables 

Vj 

species  fluid  velocity 

V 

ion/shock  velocity  ratio 

Via 

ion-acoustic  velocity 

Vi 

species  thermal  velocity 

W 

thermal  energy 

X 

spatial  coordinate 

Vi 

species  shock-fixed  velocity 

ZL 

electron  impact  ionization  rate 

Znet 

net  ionization  rate 

Zrecomb 

ion-electron  recombination  rate 

a 

fractional  ionization 

Oiia 

fractional  ionization  of  ion-acoustic  conditions 

OiS 

local  fractional  ionization  at  shock  front 

Pe 

recombination  coefficient 

7; 

species  ratio  of  specific  heats 

8 

electron  thermal  region  width  parameter 

A  Te 

electron  temperature  rise 

to 

permittivity  of  free  space 

c 

net  charge  density 
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A  neutral-neutral  mean  free  path 

A jk  inter  species  mean  free  path 

p,  Avramenko  density  parameter 

Vjk  inter-species  collision  frequency 

£  shock-centered  spatial  coordinate 

£  shock  width  parameter 

£0  charged  precursor  width 

p-  species  mass  density 

ae  electron  impact  ionization  cross  section 

Gjk  inter-species  collision  cross  section 

Tjk  inter-species  collision  time 

Vave  average  molecular  velocity 

Vk  equivalent  ionization  potential  velocity 

0  electric  potential 

uj  ion-acoustic  frequency 

u)j  species  plasma  frequency 
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Appendix  B.  Time- dependent  three-fluid  approximation 

The  preferred  method  of  analysis  of  shock  propagation  in  weakly  ionized  gases  is 
by  way  of  a  three  fluid  approach,  in  which  the  continuity,  momentum,  and  energy 
equations  of  neutrals,  ions,  and  electrons  are  numerically  solved;  the  development 
of  this  method  is  outlined  herein.  The  present  research  effort  attempted  to  lift 
the  two-fluid  approximation;  however,  the  difference  in  computational  time  steps  for 
electrons  and  heavy  particles  proved  to  be  prohibitively  large  to  be  accomplished 
within  the  time  and  computing  limits  of  this  research  effort.  In  the  Hilbun’s  two- 
fluid  solution,  ion  and  neutral  velocities  were  normalized  by  the  ion-neutral  thermal 
velocity  (22:202)  and  the  computational  time  step  was  on  the  order  of 

At2  <  Ncfl —  (87) 

Via 

where  NCfl  is  the  Courant-Friedrichs-Levy  number,  which  is  generally  0.4  and  is 
used  to  ensure  stability  (22:219).  The  three-fluid  approach  utilized  Hilbun’s  adap¬ 
tation  of  the  FCT  algorithm  of  Toth  and  Odstrcil  (38).  Therefore,  the  three-fluid 
computational  time  step  is  of  the  same  form,  but  incorporates  the  electron  thermal 
velocity  and  is  on  the  order  of 

fa  10-3At2  (88) 

Therefore,  what  had  once  taken  an  hour  to  execute  would  take  several  weeks.  It  is 
possible,  however,  to  reduce  the  time  requirements  by  a  factor  of  three.  Since  the 
motion  of  the  heavy  particles  varies  slowly  relative  to  the  electrons,  heavy  particle 
variables  only  need  to  be  updated  for  every  thousand  electron  computational  steps. 
Unfortunately,  this  still  yields  execution  times  on  the  order  of  weeks. 
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As  mentioned  above,  the  three-fluid  approach  adopts  Hilbun  s  approach  of 
solving  the  conservation  equations  for  three  species. 

^  +  ^  =  +S2  +  S3  +  SA  (89) 

C/t  C/X 

For  three  fluids,  the  vector  of  conserved  variables,  U,  becomes 


U  = 


Pn 

PnVn 

p  (ly2  +  , 1  kB—) 

rn\2v7i  '  1— 7  m  / 

Pi 

PiVi 

PiiW  +  Th-^) 

Pe 

PeVe 

0  (±v2  +  1 

re\ 2  ve  ~  l-7  me  / 


The  flux  vector,  F,  is  given  by 


F  = 


PvYn 


1-7 

r. 

% 

m  > 

r2 

i 

1 

1-7 

4 

m  ' 

n 

e 

1 

hsls.\ 

(90) 


(91) 
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The  source  terms  for  the  Euler  equations  and  species  coupling  follow  those  of  the 
two-fluid  approximation: 


0 

0 

Pn 

Pni  Pne 

Qni  “1“  Qne 

0 

0 

Pi 

,  S2  = 

1 

f3 

1 

N. 

PV5 

Qni  Qei 

0 

0 

Pe 

—Pne  +  Pei 

PeVe  _ 

Qne  “b  Qei 

(92) 


where  Pjk  represents  the  momentum  gained  by  species  j  at  the  expense  of  species 
k,  and  Qjk  represents  a  similar  transfer  of  energy  between  the  two  species.  The 
momentum  and  energy  coupling  terms  are  given  by  Jaffrin  (23:611): 


Pni  =  |nnnicrin(Vri  -  Vn)y^-kB(Tn  +  T») 


(93) 


D  _  8  /t/  T/  X  l2mekBTe 

Pne  —  ^Tln7ie<7eri(Ve  Ki)  y  ^ 


(94) 


8  trN  /2meA :BTe 

Pei  =  - nineoei(Vi  -  Ve)y - - - 


(95) 


Qni  =  2nnniain\  f-^—kB(Tn  +  Ti){kB(Tn  -  Pi)  +  om(^  ~  ^)(^  +  ^n)}  (96) 

V  7rm  <1 
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Q„  =  8W, -  T„)  +  |(V.  -  V.)(K,  +  K) }  (97) 


Qei  =  Snineae 


1 2 mekBTe  ,  kB 


mP 


4^(7)  _  T.)  +  i(V5  -  K)(V)  +  -f  K)}  (98) 

7T  TO  O  m 


The  electric  field  source  term  is  given  by 


S3  =  E 


0 

0 

0 

0 

ern 

eriiVi 

0 

-ene 

-eneVe 


(99) 


where  E  is  the  local  electric  field  and  e  is  the  elementary  charge.  Finally,  the 
ionization  source  term  is  given  by 


54  = 


Znet'm 

Znet'm 

ZnetmVi 

Z^m{\V?  + 

Z-net^e 


ZnetmeYe 

Zmtme(lV?  +  ^) 


(100) 
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In  order  to  prevent  the  plasma  from  decaying  in  the  case  of  the  ternary  fluid  ap¬ 
proximation  of  the  Riemann  problem,  it  is  either  necessary  to  apply  the  restriction 
of  Znet  >  0  or  apply  an  ambient  electric  field  to  maintain  the  plasma.  Hilbun’s 
two-fluid  code  is  easily  expanded  to  include  the  electron  equations.  Nine  variables- 
density,  velocity,  and  temperature  for  each  of  the  components-are  solved  using  these 
nine  equations.  Ideal  behavior  of  the  three  components  is  assumed  to  determine 
pressures. 

The  flow  variables  were  non-dimensionalized  by  the  following  normalizing  pa¬ 
rameters:  the  ambient  upstream  (Region  1)  neutral  density,  pnQ)  the  ambient  up¬ 
stream  neutral  and  electron  temperatures,  Tn o  and  Te o,  respectively,  the  neutral  and 
electron  thermal  velocities  at  their  ambient  temperatures,  vn  and  ve,  respectively, 
and  the  length  of  the  shock  tube,  L.  Using  these  parameters,  the  normalized  vari¬ 
ables  are: 


x'  = 


x 

L’ 


,  rne  pe 
Pe  ^ 

m  PnO 


T 

rpt  _  in 
^  n  rri  > 

inO 


(101) 


In  the  three  fluid  solution,  the  steady-state  electric  field  approximation  (Equa¬ 
tion  48)  is  no  longer  valid.  Therefore,  Poisson’s  equation,  Equation  2,  must  be 
solved  directly  at  each  point  in  space  in  order  to  obtain  the  field.  Recall  the  one 
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dimensional  form  of  Poisson’s  equation: 


dE  _  e(nj  -  ne) 
dx  6o 

Since  the  problem  is  one-dimensional,  the  net  charge  density  at  each  point  can  be 
treated  as  an  infinite  sheet  of  charge  (infinite  in  the  directions  normal  to  the  prob¬ 
lem’s  x-direction).  With  this,  the  electric  field  contributions  form  each  point  are 
independent  of  distance.  The  field  evaluates  to 
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ionization  and  ionization  mitigated  by  ion-electron  recombination.  Introduction  of  a  spatial  variation  in  electron  temperature 
resulted  in  a  broadening  and  strengthening  of  the  electric  field  associated  with  the  electronic  double  layer  (EDL)  at  the  shock 
front.  Results  of  unrestricted  ionization  were  a  broadening  and  strengthening  of  the  electric  field  associated  with  the  EDL,  an 
acceleration  of  the  neutral  shock  front,  and  the  development  of  a  neutral  precursor  ahead  of  the  shock.  Ion-electron 
recombination  was  seen  to  reduce  these  effects. 
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